﻿#pragma once

#ifndef __MELTING_MATH_H__
#define __MELTING_MATH_H__

//--------------------------------------------------------------------------------------------------------------
// platform 
#ifdef WIN32	
	
	#include <math.h>
	
	#define WIN32_LEAN_AND_MEAN	// do not include unnecessary headers
	#include <windows.h>

	//#define _XM_NO_INTRINSICS_	// by defining '_XM_NO_INTRINSICS_', we allow the same Win32 & XBOX360 code use of the XMVECTOR,etc.
	//#include <xnamath.h>		// NB : this temporary !

#endif

//#ifdef _XBOX
//	#include <xtl.h>
//	#include <xboxmath.h>
//#endif

#include <xmmintrin.h>
#include <algorithm>
#include <vector>
//-------------------------------------------------------------------------------------------------------------- 
#ifdef _DEBUG
#else
#endif

#ifndef LOG
#	define LOG		printf
#endif

#ifndef MB_ASSERT
#	define MB_ASSERT( expr, ... )					if( !(expr) ) {__debugbreak();}
#endif 

#ifndef UTEST( expression )
#	define UTEST( expression ) 		if( !(expression ) )			{	LOG("- TEST FAILED : "#expression"\n" ); bTestSuccess = false; } 
#endif 
#ifndef UTEST_START					
#	define UTEST_START				bool bTestSuccess = true;
#endif 
#ifndef UTEST_RESULT
#	define UTEST_RESULT				bTestSuccess
#endif 
 

#define ACTIVATE_ASSERT_ON_UNIT_LENGTH_VECTOR 1

#if ACTIVATE_ASSERT_ON_UNIT_LENGTH_VECTOR
#	define 	ASSERT_UNIT_LENGTH_VECTOR(_vVector) 	MB_ASSERT( Math::Abs( REAL( 1.0f) - _vVector.GetLength() ) < REAL( 0.01f ) , "Vector "#_vVector " is not unit length\n" ); 
#else
# 	define 	ASSERT_UNIT_LENGTH_VECTOR(_vVector)		
#endif	
//--------------------------------------------------------------------------------------------------------------
 
typedef		float			REAL;
typedef		REAL			Radian;
typedef		REAL			Degree;

typedef		float	FLOAT; 

struct Float2
{
	union
	{
		struct
		{ 
			FLOAT x;
			FLOAT y;
		};
		FLOAT a[ 2 ];
	};

	Float2( FLOAT _x, FLOAT _y ) 
	: x( _x )
	, y( _y )
	{}
};

typedef		unsigned char	UINT8;
typedef		unsigned short	UINT16;
typedef		unsigned int	UINT32;	 
 
typedef		signed char		INT8;
typedef		signed short	INT16;
typedef		signed int		INT32;

#define		VectorRight()	Vec3(1,0,0)
#define		VectorLeft()	Vec3(-1,0,0)
#define		VectorUp()		Vec3(0,1,0)
#define		VectorDown()	Vec3(0,-1,0)
#define		VectorFront()	Vec3(0,0,1)
#define		VectorBack()	Vec3(0,0,-1)
	
//--------------------------------------------------------------------------------------------------------------


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	Declarations :	UTILS - Utility functions for math ,3D, and other related stuff
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting 
{
	namespace Math
	{  
		static const REAL PI	= 3.141592654f;			///< PI

		/// Are the two given REAL equal ? with respect to the given epsilon 
		inline static bool					Equal( REAL _realA, REAL _realB, REAL _epsilon = 0.001f ); 

		/// Get the inverse of the given REAL
		inline static REAL					Inverse	( REAL _real );
		
		/// Swap the 2 values in the given parameters
		inline static void 					Swap( REAL & _rVal1, REAL & _rVal2 );			

		/// Clamp the given value between min and max
		inline static void					Clamp( REAL & _rfOutClamped, REAL _fMin, REAL _fMax );

		/// Get the angle in Radian matching the one given in Degree
		inline static Radian				DegreeToRadian( Degree _fAngle );
		/// Get the angle in Degree matching the one given in Radian
		inline static Degree				RadianToDegree( Radian _fAngle );


		/// Get the Sin of the given angle 
		inline static REAL 					Sin	( Radian _fAngle );
		/// Get the ArcSin of the given value, in the range –π/2 to π/2 radians
		inline static Radian 				ASin( REAL _fValue );
		/// Get the Cos of the given angle
		inline static REAL 					Cos	( Radian _fAngle );
		/// Get the ArcCos of the given value, in the range 0 to π radians
		inline static Radian 				ACos( REAL _fValue );
		/// Get the Tan of the given angle 
		inline static REAL 					Tan	( Radian _fAngle );
		/// Get the ArcTan of the given value, in the range of –π/2 to π/2 radians 
		inline static Radian 				ATan( REAL _fValue );
		/// Get the ArcTan2 of _fValueY / _fValueX (taking the signs into consideration), in the range –π to π radians
		inline static Radian 				ATan2( REAL _fValueY, REAL _fValueX );
		 
		/// Get the Sqr of the given REAL
		inline static REAL					Square( REAL _real );			
		/// Get the Sqrt of the given REAL
		inline static REAL					SquareRoot( REAL _real );	
		/// Get the 1.f / Sqrt of the given REAL
		inline static REAL					InverseSquareRoot( REAL _real );
		/// Get absolute value
		inline static REAL					Abs( REAL _real );

		/// Calculate the smallest multiple that is equal to or greater than the specified value.
		inline static UINT32				NextMultiple( UINT32 _ulValue, UINT32 _ulMultiple );

		/// Unit testing every functions in this section
		inline bool							UtilsUnitTest( void );

	}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	Declarations :	VECTOR - multi dimensional crap --- to be re-coded !
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting
{  
	namespace Math
	{
		namespace
		{
			template< class Component >
			class AffectationPolicy
			{
			public:
				static void Affect( Component& c, REAL fValue ) { c = Component( fValue.Round() ); }
			};
			
			template<>
			class AffectationPolicy<REAL>
			{
			public:
				static void Affect( REAL& c, REAL fValue ) { c = fValue; }
			};
		}
		
		//////////////////////////////////////////////////////////
		//  class Vector
		/// Mutidimensional vector with common operations on it (length, dot product, etc.).
		//////////////////////////////////////////////////////////
		template< class Component, UINT32 Dim >
		class Vector
		{	
		protected:
			Component				mtVector[Dim];
		
		public:
			typedef Component ComponentType;
			enum { Dimension = Dim };
		
			Vector( void );
			///	With this Ctor, val is used to initialize each component of the Vector
			Vector( const Component* aVec );		
			/// Every component of the Vector is set to 'value'		
			explicit Vector( const Component& cmp );
			Vector( const Vector<Component, Dim >& vecOther );
			~Vector( void )									{ }
			
			inline Vector<Component, Dim>&							operator=( const Vector<Component, Dim>& vecOther );
			
			inline Vector<Component, Dim> 							operator-( void ) const	{ return Vector<Component, Dim>( Component(0) ) -= *this; }
			
			inline Vector<Component, Dim>  							operator+( const Vector<Component, Dim>& vecOther ) const	{ return Vector<Component, Dim>( *this ) += vecOther; }
			inline Vector<Component, Dim> 		 					operator-( const Vector<Component, Dim>& vecOther ) const	{ return Vector<Component, Dim>( *this ) -= vecOther; }
			inline Vector<Component, Dim>  							operator*( const Vector<Component, Dim>& vecOther ) const	{ return Vector<Component, Dim>( *this ) *= vecOther; }
			inline Vector<Component, Dim>  							operator/( const Vector<Component, Dim>& vecOther ) const	{ return Vector<Component, Dim>( *this ) /= vecOther; }
			inline Vector<Component, Dim>  							operator*( const Component& cmp ) const						{ return Vector<Component, Dim>( *this ) *= Vector<Component, Dim>( cmp ); }
			
			inline Vector<Component, Dim>&							operator+=( const Vector<Component, Dim>& vecOther );
			inline Vector<Component, Dim>&							operator-=( const Vector<Component, Dim>& vecOther );
			inline Vector<Component, Dim>&							operator*=( const Vector<Component, Dim>& vecOther );
			inline Vector<Component, Dim>&							operator*=( const REAL& fxValue );
			inline Vector<Component, Dim>&							operator/=( const Vector<Component, Dim>& vecOther );
			inline Vector<Component, Dim>&							operator/=( const REAL& fxValue );
			
			inline bool										operator==( const Vector& vecOther ) const;
			inline bool										operator!=( const Vector& vecOther ) const					{ return !( *this == vecOther ); }
			
			static inline void								LERP( REAL _fTime, const Vector<Component, Dim>& rvA, const Vector<Component, Dim>& rvB, Vector<Component, Dim>& rvOutResult );
			static inline void 								CrossP( const Vector<Component, Dim>& rvA, const Vector<Component, Dim>& rvB, Vector<Component, Dim>& rvOutResult );
			inline Vector<Component, Dim >	 				CrossP( const Vector<Component, Dim>& vecOther ) const;
			inline REAL										DotP( const Vector<Component, Dim>& vecOther ) const;
			inline REAL										GetLength( void ) const;
			inline REAL										GetSquaredLength( void ) const;
			inline Vector<Component, Dim>&					Normalize( void );
			inline Vector<Component, Dim>&					NormalizeInteger( void );
			
			inline void										SetComponent( UINT32 ulIndex, const Component& tVal );
			inline Component								GetComponent( UINT32 ulIndex ) const;
			inline Component&								operator[]( UINT32 ulIndex );
			inline const Component&							operator[]( UINT32 ulIndex ) const;


		// TODO : put the func below in Math::
			/// Return wether this vector is within the given range(min,max)
			inline bool										IsWithinRange( const Vector<Component, Dim>& rvMin, const Vector<Component, Dim>& rvMax );
			
			/// Return in the reference parameter 'rvOutResult' a vector composed foreach index 'i' of the min(rvA[i], rvB[i])
			inline static void								GetVectorMin( const Vector<Component, Dim>& rvA, const Vector<Component, Dim>& rvB, Vector<Component, Dim>& rvOutResult );

			/// Return in the reference parameter 'rvOutResult' a vector composed foreach index 'i' of the max(rvA[i], rvB[i])
			inline static void								GetVectorMax( const Vector<Component, Dim>& rvA, const Vector<Component, Dim>& rvB, Vector<Component, Dim>& rvOutResult );
			
			/// Return in the reference parameter 'rvOutResult' a vector composed foreach index 'i' of the 1/rvA[i]
			inline static void								GetVectorInverse( const Vector<Component, Dim>& rvA, Vector<Component, Dim>& rvOutResult );
				
	#ifndef __FINALROM 
			inline static bool								UnitTest( void );
			inline void										Print( void ) const;
	#endif
		};
		
		
		
		//////////////////////////////////////////////////////////
		//  class Vector2
		/// Helper class that provides custom ctors, operations, etc., for bidimensional vectors.
		//////////////////////////////////////////////////////////
		template< class Component >
		class Vector2 : public Vector< Component, 2 >
		{
		public: 
		
			explicit Vector2( const Component& cmp )			: Vector< Component, 2 >( cmp )			{}
			Vector2( const Component& x, const Component& y )	: Vector< Component, 2 >()				{ SetX( x ); SetY( y ); }
			Vector2( const Vector< Component, 2 >& vectOther )	: Vector< Component, 2 >( vectOther )	{}
			Vector2( const Vector2< Component >& vectOther )	: Vector< Component, 2 >( vectOther )	{}
			Vector2( void )										: Vector< Component, 2 >()				{}
		
			inline Component			GetX() const															{ return (*this)[ 0 ]; }
			inline Component			GetY() const															{ return (*this)[ 1 ]; }
			inline void					SetX( const Component& x )												{ (*this)[ 0 ] = x; }
			inline void					SetY( const Component& y )												{ (*this)[ 1 ] = y; }
			inline void					SetXY( const Component& x, const Component& y )							{ SetX( x ); SetY( y ); }
			
			Radian						GetAngleWith( const Vector2< Component >& vOther ) const;
			Vector2<Component>&			Rotate( Radian fAngle );
			
			REAL						GetSquaredDistanceToSegment( const Vector2< Component >& vA, const Vector2< Component >& vB ) const;
			REAL						GetDistanceToSegment( const Vector2< Component >& vA, const Vector2< Component >& vB ) const;		
			
			static const Vector2< Component > ZERO;
			static const Vector2< Component > ONE;
			
			inline static bool								UnitTest( void );

		private:
		};
		
		template< typename Component > 
		const Vector2< Component > Vector2< Component >::ZERO ( Component(0), Component(0));
		
		template< typename Component > 
		const Vector2< Component > Vector2< Component >::ONE ( Component(1), Component(1));
		
		typedef Vector2<REAL>	Vec2;
		typedef Vector2<REAL>	Vec2REAL;
		typedef Vector2<INT32>	Vec2INT32;
		typedef Vector2<INT16>	Vec2INT16;
		typedef Vector2<UINT32>	Vec2UINT32;
		typedef Vector2<UINT16>	Vec2UINT16;
		typedef Vector2<INT8>	Vec2INT8;
		typedef Vector2<UINT8>	Vec2UINT8;
		
		
		//////////////////////////////////////////////////////////
		//  class Vector2 implemntation
		//////////////////////////////////////////////////////////
		template< class Component, UINT32 Dim >
		Vector<Component, Dim >::Vector( void )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] = (Component)0;
			}
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim >::Vector( const Component& val )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] = val;
			}
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim >::Vector( const Component* aVec )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] = aVec[i];
			}
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim >::Vector( const Vector<Component, Dim>& vecOther )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] = vecOther.mtVector[i];
			}
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim >& Vector<Component, Dim >::operator=( const Vector<Component, Dim>& vecOther )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] = vecOther.mtVector[i];
			}
			
			return *this;
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim >& Vector<Component, Dim>::operator+=( const Vector<Component, Dim>& vecOther )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] += vecOther.mtVector[i];
			}
			
			return *this;
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim >& Vector<Component, Dim>::operator-=( const Vector<Component, Dim>& vecOther )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] -= vecOther.mtVector[i];
			}
			
			return *this;
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim >& Vector<Component, Dim>::operator*=( const Vector<Component, Dim>& vecOther )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] *= vecOther.mtVector[i];
			}
			
			return *this;
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim>& Vector<Component, Dim>::operator/=( const Vector<Component, Dim>& vecOther )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] /= vecOther.mtVector[i];
			}
			
			return *this;
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim >& Vector<Component, Dim>::operator*=( const REAL& fxValue )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] *= fxValue;
			}
			
			return *this;
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim >& Vector<Component, Dim>::operator/=( const REAL& fxValue )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				mtVector[i] /= fxValue;
			}
			
			return *this;
		}
		
		template< class Component, UINT32 Dim >
		bool Vector<Component, Dim >::operator==( const Vector<Component, Dim>& vecOther ) const
		{
			REAL fEpsilon = 0.00001f;

			for( UINT32 i = 0; i < Dim; i++ )
			{
				if( Math::Abs( mtVector[i] - vecOther.mtVector[i]) >= fEpsilon )
				{
					return false;
				}
			}
			
			return true;
		}
		
		template< class Component, UINT32 Dim >
		void Vector<Component, Dim>::LERP( REAL _fTime, const Vector<Component, Dim>& rvA, const Vector<Component, Dim>& rvB, Vector<Component, Dim>& _rvOutResult ) 
		{
			_rvOutResult = rvA + (rvB - rvA) * _fTime;
		}
		
		template<>
		void Vector<REAL, 3>::CrossP( const Vector< REAL, 3 >& rvA, const Vector< REAL, 3 >& rvB, Vector< REAL, 3 >& rvOutResult ) 
		{
			rvOutResult[0] = rvA.mtVector[1] * rvB.mtVector[2] - rvB.mtVector[1] * rvA.mtVector[2];
			rvOutResult[1] = rvA.mtVector[2] * rvB.mtVector[0] - rvB.mtVector[2] * rvA.mtVector[0];
			rvOutResult[2] = rvA.mtVector[0] * rvB.mtVector[1] - rvB.mtVector[0] * rvA.mtVector[1];
		}
		
		// NDAL : sorry I'll simply make it blow for n  (!=3) dimensions
		template< class Component, UINT32 Dim >
		void Vector<Component, Dim>::CrossP( const Vector<Component, Dim>& rvA, const Vector<Component, Dim>& rvB, Vector<Component, Dim>& rvOutResult ) 
		{
			MB_ASSERT(false, "ARGL!! Who needs a CrossProduct for nDimension vectors ???");
		}
		
		// Cross Product for 3D
		template<>
		Vector<REAL, 3> Vector<REAL, 3>::CrossP( const Vector<REAL, 3 >& vecOther ) const
		{
			REAL fxVec[3] = 
			{
				mtVector[1] * vecOther.mtVector[2] - vecOther.mtVector[1] * mtVector[2],
				mtVector[2] * vecOther.mtVector[0] - vecOther.mtVector[2] * mtVector[0],
				mtVector[0] * vecOther.mtVector[1] - vecOther.mtVector[0] * mtVector[1]
			};
			
			return Vector<REAL, 3>( fxVec );
		}

		// for nD, its a pain in the ass to compute a matrix determinant, and it's not really useful
		// so returns default vector
		// ...NDAL : but we'll make it blow with an assert just in case someone might be tempted to use it...
		template< class Component, UINT32 Dim >
		Vector<Component, Dim> Vector<Component, Dim>::CrossP( const Vector<Component, Dim>& vecOther )  const
		{
			MB_ASSERT(false, "ARGL!! Who needs a CrossProduct for n (!=3) dimensions vectors ???");
			
			return Vector<Component, Dim>();
		}
		
		template< class Component, UINT32 Dim >
		REAL Vector<Component, Dim >::DotP( const Vector<Component, Dim>& vecOther ) const
		{
			REAL fResult( 0 );
			
			for( UINT32 i = 0; i < Dim; i++ )
			{
				fResult += REAL( mtVector[i] * vecOther.mtVector[i] );
			}
			
			return fResult;
		}
		
		template< class Component, UINT32 Dim >
		REAL Vector<Component, Dim >::GetLength( void ) const
		{  
			REAL fResult( 0 );
			
			for( UINT32 i = 0; i < Dim; i++ )
			{
				REAL fValue(mtVector[i]);
				fResult += REAL( fValue * fValue );
			}
			
			return Math::SquareRoot( fResult );
		}
		
		template< class Component, UINT32 Dim >
		REAL Vector<Component, Dim >::GetSquaredLength( void ) const
		{
			REAL fResult( 0 );
			
			for( UINT32 i = 0; i < Dim; i++ )
			{
				fResult += REAL( mtVector[i] * mtVector[i] );
			}
			
			return REAL( fResult );
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim>& Vector<Component, Dim >::Normalize( void )
		{
			REAL fxLength( GetLength() );
					
			if( ( 1 != fxLength ) && ( 0 != fxLength ) )
			{ 
				REAL fxInvLength = 1.f / fxLength;
				
				for( UINT32 i = 0; i < Dim; i++ )
				{
					REAL resultValue( mtVector[i] );
					resultValue *= fxInvLength;
					mtVector[i] = Component( resultValue );
				}
			} 		
		
	//		MB_ASSERT( Math::Abs( 1 - GetLength() ) < 0.01f, "Could not normalize the Vector T_T" );		 

			return *this;
		}
		
		template< class Component, UINT32 Dim >
		Vector<Component, Dim>& Vector<Component, Dim >::NormalizeInteger( void )
		{
			REAL fxLength( GetLength() );
			
			if( ( 1 == fxLength ) || ( 0 == fxLength ) )
			{ 
				return *this;;
			}
			
			// old way		
	#if 0		
			for( UINT32 i = 0; i < Dim; i++ )
			{
				REAL resultValue( mtVector[i] );
				
				resultValue /= fxLength;
				
				mtVector[i] = Component( resultValue.Round() );
				
			}
	#endif		
			REAL fxInvLength = fxLength.Inverse();
			
			for( UINT32 i = 0; i < Dim; i++ )
			{
				REAL resultValue( mtVector[i] );
				resultValue *= fxInvLength;
				mtVector[i] = Component( resultValue.Round() );
			}
			
			return *this;
		}
		
		template< class Component, UINT32 Dim >
		inline void Vector<Component, Dim >::SetComponent( UINT32 ulIndex, const Component& tVal )
		{
			mtVector[ulIndex] = tVal;
		}
		
		template< class Component, UINT32 Dim  >
		inline Component Vector<Component, Dim >::GetComponent( UINT32 ulIndex ) const
		{
			return Component(mtVector[ulIndex]);
		}
			
	#ifndef __FINALROM 
		
		template< class Component, UINT32 Dim >
		inline void	Vector<Component, Dim>::Print( void ) const
		{
			LOG("(%f,%f,%f)\n", mtVector[0], mtVector[1], mtVector[2] );

			for( UINT32 i = 0; i < Dim; i++ )
			{//mtVector[i].Print();
			}
		}
	#endif
		
		template< class Component, UINT32 Dim >
		inline Component& Vector<Component, Dim>::operator[]( UINT32 ulIndex )
		{
			return mtVector[ ulIndex ];
		}
		
		template< class Component, UINT32 Dim >
		inline const Component& Vector<Component, Dim>::operator[]( UINT32 ulIndex ) const
		{
			return mtVector[ ulIndex ];		
			
		}


		// Return wether this vector is within the given range(min,max)
		template< class Component, UINT32 Dim >
		bool Vector<Component, Dim>::IsWithinRange( const Vector<Component, Dim>& rvMin, const Vector<Component, Dim>& rvMax )
		{
			bool bIsWithingRange = true;
			
			for( UINT32 i = 0; i < Dim; i++ )
			{
				if( rvMin.mtVector[ i ] > mtVector[ i ] || mtVector[ i ] > rvMax.mtVector[ i ] )
				{
					bIsWithingRange = false;
					break;
				}
			}
			
			return bIsWithingRange;
		}
		
		// Return in the reference parameter 'rvOutResult' a vector containing foreach component index 'i' the min(rvA[i], rvB[i])
		template< class Component, UINT32 Dim >
		void Vector<Component, Dim>::GetVectorMin( const Vector<Component, Dim>& rvA, const Vector<Component, Dim>& rvB, Vector<Component, Dim>& rvOutResult )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				rvOutResult [ i ] = min( rvA[ i ], rvB[ i ]);
			}
		}

		// Return in the reference parameter 'rvOutResult' a vector containing foreach component index 'i' the max(rvA[i], rvB[i])
		template< class Component, UINT32 Dim >
		void Vector<Component, Dim>::GetVectorMax( const Vector<Component, Dim>& rvA, const Vector<Component, Dim>& rvB, Vector<Component, Dim>& rvOutResult )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				rvOutResult [ i ] = max( rvA[ i ], rvB[ i ]);
			}
		}
		
		// Return in the reference parameter 'rvOutResult' a vector containing foreach component index 'i' the max(rvA[i], rvB[i])
		template< class Component, UINT32 Dim >
		void Vector<Component, Dim>::GetVectorInverse( const Vector<Component, Dim>& rvA, Vector<Component, Dim>& rvOutResult )
		{
			for( UINT32 i = 0; i < Dim; i++ )
			{
				rvOutResult [ i ] =  rvA[ i ];
			}
		}	 
		
		template< class Component >
		inline Vector<Component, 2> MakeVector2D( const Component& tComp1, const Component& tComp2 )
		{
			Vector< Component, 2 > vecResult;
			vecResult[0] = tComp1;
			vecResult[1] = tComp2;

			return vecResult;
		}
		
		
		
		template< class Component >
		Radian Vector2< Component >::GetAngleWith( const Vector2< Component >& vOther ) const
		{
			REAL	fxDotProduct	= DotP( vOther );
			REAL	fxLength		= GetLength() * vOther.GetLength();
			Radian	radAngle		= 0;
			REAL	fDotPerLength	= fxDotProduct / fxLength;

			if( Math::Abs( fDotPerLength ) > 1.f )
			{
				radAngle = Radian( 0.f );
			}
			else
			{
				radAngle		= (Radian) Math::ACos( fDotPerLength );
				bool bPositive	= ( GetX() * vOther.GetY() ) >= ( GetY() * vOther.GetX() );
				radAngle		= ( bPositive ? radAngle : -radAngle );
			}
			
			return radAngle;
		}
		
		template< class Component >
		Vector2< Component >& Vector2< Component >::Rotate( Radian fAngle )
		{
			REAL fcos( Math::Cos( fAngle ) );
			REAL fsin( Math::Sin( fAngle ) );

			REAL fX( GetX() );
			REAL fY( GetY() );
			
			REAL fXNew = REAL( fX * fcos - fY * fsin );
			REAL fYNew = REAL( fX * fsin + fY * fcos );
			
			AffectationPolicy<Component>::Affect( mtVector[ 0 ], fXNew );
			AffectationPolicy<Component>::Affect( mtVector[ 1 ], fYNew );
			
			return *this;
		}
		 
		template< class Component >
		REAL Vector2< Component >::GetDistanceToSegment( const Vector2< Component >& vA, const Vector2< Component >& vB ) const
		{
			return Math::SquareRoot( GetSquaredDistanceToSegment( vA, vB ) );
		}
		
		template< class Component >
		REAL Vector2< Component >::GetSquaredDistanceToSegment( const Vector2< Component >& vA, const Vector2< Component >& vB ) const
		{
			if( vA.GetX() == vB.GetX() && vA.GetY() == vB.GetY() )
			{
				return ( *this - vA ).GetSquaredLength();
			}
			
			Component cXDist = vB.GetX() - vA.GetX();
			Component cYDist = vB.GetY() - vA.GetY();
			
    		REAL fxRatio = REAL( ( GetX() - vA.GetX() ) * cXDist + ( GetY() - vA.GetY() ) * cYDist ) /
    			REAL( cXDist * cXDist + cYDist * cYDist );
	    		
    		if( fxRatio < REAL( 0.f ) )
    		{
    			return ( *this - vA ).GetSquaredLength();
    		}
    		else if( fxRatio > 1 )
    		{
    			return ( *this - vB ).GetSquaredLength();
    		}
    		else
    		{
    			Vector2<REAL> vNew( ( 1 - fxRatio ) * REAL( vA.GetX() ) + fxRatio * REAL( vB.GetX() ),
    				( 1 - fxRatio ) * REAL( vA.GetY() ) + fxRatio * REAL( vB.GetY() ) );
	    			
    			REAL fxXMe = REAL( GetX() );
    			REAL fxYMe = REAL( GetY() );
	    		
    			Vector2<REAL> vResult;
    			vResult.SetX( fxXMe - vNew.GetX() );
    			vResult.SetY( fxYMe - vNew.GetY() );    				
    			return vResult.GetSquaredLength();
    		}
		}
		
		///////////////////////////////////////////////////////////////////////////////////////////////
		/// Helper class that provide custom ctors, operations, etc., for tridimensional vectors.
		///////////////////////////////////////////////////////////////////////////////////////////////
		template< class Component >
		class Vector3 : public Vector< Component, 3 >
		{
		public:
			explicit Vector3( const Component& cmp )	: Vector< Component, 3 >( cmp )	{}
			Vector3( void )								: Vector< Component, 3>()		{}
			Vector3( const Component& x, const Component& y, const Component& z )		{ mtVector[0] = x; mtVector[1] = y; mtVector[2] = z; }
			Vector3( const Vector< Component, 3 >& other )								{ mtVector[0] = other[ 0 ]; mtVector[1] = other[ 1 ]; mtVector[2] = other[ 2 ]; }
		
			inline Component	GetX() const										{ return (*this)[ 0 ]; }
			inline Component	GetY() const										{ return (*this)[ 1 ]; }
			inline Component	GetZ() const										{ return (*this)[ 2 ]; }
			
			inline void			SetX( const Component& x )							{ (*this)[ 0 ] = x; }
			inline void			SetY( const Component& y )							{ (*this)[ 1 ] = y; }
			inline void			SetZ( const Component& z )							{ (*this)[ 2 ] = z; }
			inline void			SetXYZ( const Component& x, const Component& y, const Component& z )		{ SetX( x ); SetY( y ); SetZ( z ); }
			Vector3<Component>&	RotateXZ( REAL fAngle );
		
			Radian				GetAngleWith( const Vector3< Component >& vOther ) const;
			
			static const Vector3< Component > ZERO ;
			static const Vector3< Component > ONE ;
			
		private:
		};
		
		template< typename Component > 
		const Vector3< Component >  Vector3< Component >::ZERO ( Component(0), Component(0), Component(0) );
		template< typename Component > 
		const Vector3< Component >  Vector3< Component >::ONE ( Component(1), Component(1), Component(1) );
		
		typedef Vector3<REAL>	Vec3;
		typedef Vector3<REAL>	Vec3REAL;
		typedef Vector3<INT32>	Vec3INT32;
		typedef Vector3<INT16>	Vec3INT16;
		typedef Vector3<UINT32>	Vec3UINT32;
		typedef Vector3<UINT16>	Vec3UINT16;
		typedef Vector3<INT8>	Vec3INT8;
		typedef Vector3<UINT8>	Vec3UINT8;
		
		template< class Component >
		Vector3< Component >& Vector3< Component >::RotateXZ( REAL fAngle )
		{
			
			REAL fX( GetX() );
			REAL fZ( GetZ() );
			
			Vector2<REAL> tempVect( fX , fZ );
			
			tempVect.Rotate( fAngle );
			
			AffectationPolicy<Component>::Affect( mtVector[ 0 ], tempVect.GetX() );
			AffectationPolicy<Component>::Affect( mtVector[ 2 ], tempVect.GetY() );
			
			return *this;
		}
		
		template< class Component >
		Radian Vector3< Component >::GetAngleWith( const Vector3< Component >& vOther ) const
		{
			REAL fScalar( DotP( vOther ) );
			REAL fNorms( GetLength() * vOther.GetLength() );
			REAL fCos( fScalar / fNorms );

			return Math::Abs( fCos ) > 1.0f ? 0.f : Math::ACos( fCos );
		}
	}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	Declarations :	MATRIX --- to be re-coded
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting
{
	namespace Math
	{

	#define		ROUND_UP(_value,_align)		( _value + ( _align - 1 ) & (~( _align - 1) )

	// for MSVC align
	#define		ALIGN(_stuff,_align)		__declspec(align(_align))	_stuff
	// for GCC align
	//#define	ALIGN(_stuff,_align)		_stuff		__attribute__ ((aligned (_align))).

	#define		ALIGN_16(_stuff)		ALIGN(_stuff,16)

 
		class Quaternion;
		class SRT; 

		/*
			TODO : use branch less code with _mm_cmplt_ps and _mm_blendv_ps 

			xmmDot = SSEDotProduct( ... );

			xmmZero = _mm_setzero_ps();


			xmmMask = _mm_cmplt_ps( xmmDot, xmmZero );
			xmmTrue = _mm_load_ps( & valueTrue[ ... ] );
			xmmFalse = _mm_load_ps( & valueFalse[ ... ] );

			xmmResult = _mm_blendv_ps(xmmTrue, xmmFalse, xmmMask );  
		*/

#define SIMD_PREFETCH_NOREUSE(mem)			_mm_prefetch((char*)(mem),_MM_HINT_NTA)
		

#define SIMD_LOAD_A_FLOAT_AS_000X( float1, xmm )											\
					{(xmm) = _mm_set_ps( 0.0f, 0.0f, 0.0f, (float1) );}		
#define SIMD_LOAD_A_FLOAT_AS_XXXX( float1, xmm )											\
					{(xmm) = _mm_set_ps( (float1), (float1), (float1), (float1) );}		

 // OK - but since it loads like this (z,z,y,x), it is not appropriate for every processing - e.g. summing the component is wrong!
#define SIMD_LOAD_A_FLOAT3_AS_ZZYX( float3, xmm )											\
					{(xmm) = _mm_loadu_ps( (float*)& (float3) );							\
					 (xmm) = _mm_shuffle_ps( (xmm), (xmm), _MM_SHUFFLE(2,2,1,0) );}		

#define SIMD_LOAD_A_FLOAT3_AS_0ZYX( float3, xmm )											\
					{(xmm) = _mm_set_ps( 0.0f, (float3)[2], (float3)[1], (float3)[0] );}		

#define SIMD_LOAD_A_FLOAT3_AS_1ZYX( float3, xmm )											\
					{(xmm) = _mm_set_ps( 1.0f, (float3)[2], (float3)[1], (float3)[0] );}		

#define SIMD_LOAD_A_FLOAT3	SIMD_LOAD_A_FLOAT3_AS_0ZYX

#define SIMD_LOAD_A_FLOAT4( float4, xmm )											\
					{(xmm) = _mm_set_ps( (float4)[3], (float4)[2], (float4)[1], (float4)[0] );}		

#define SIMD_LOAD_A_4FLOAT( float1X, float1Y, float1Z, float1W, xmm )				\
					{(xmm) = _mm_set_ps( float1W, float1Z, float1Y, float1X );}		


#define SIMD_STORE_A_FLOAT3( float3, xmm )													\
					{ ALIGN_16( float tmp[4] );												\
					 _mm_store_ps( tmp, (xmm) );											\
					 (float3)[0] = tmp[0]; (float3)[1] = tmp[1]; (float3)[2] = tmp[2]; }

#define SIMD_STORE_A_FLOAT4( float4, xmm )													\
					{ _mm_store_ps( float4, (xmm) );}

#define SIMD_SUM_COMPONENTS( a, r )															\
					{(r) = _mm_add_ps((a),_mm_movehl_ps((a),(a)) ); 						\
					 (r) = _mm_add_ps((r), _mm_shuffle_ps( (r), (r), _MM_SHUFFLE(1,1,1,1))); }

#define	SIMD_DOT_PRODUCT( a, b, r )							\
					{SIMDVector ab = _mm_mul_ps((a),(b));	\
					 SIMD_SUM_COMPONENTS(ab,(r));}

#define	SIMD_CROSS_PRODUCT_3( a, b, r )																											\
					{SIMDVector c = _mm_mul_ps(_mm_shuffle_ps((a),(a),_MM_SHUFFLE(0,0,2,1)), _mm_shuffle_ps((b),(b),_MM_SHUFFLE(0,1,0,2)));		\
					 SIMDVector d = _mm_mul_ps(_mm_shuffle_ps((a),(a),_MM_SHUFFLE(0,1,0,2)), _mm_shuffle_ps((b),(b),_MM_SHUFFLE(0,0,2,1)));		\
					 (r) = _mm_sub_ps(c,d);}

#define	SIMD_CROSS_PRODUCT_4( a, b, r )																											\
					{SIMDVector c = _mm_mul_ps(_mm_shuffle_ps((a),(a),_MM_SHUFFLE(0,3,2,1)), _mm_shuffle_ps((b),(b),_MM_SHUFFLE(1,0,3,2)));		\
					 SIMDVector d = _mm_mul_ps(_mm_shuffle_ps((a),(a),_MM_SHUFFLE(1,0,3,2)), _mm_shuffle_ps((b),(b),_MM_SHUFFLE(0,3,2,1)));		\
					 (r) = _mm_sub_ps(c,d);}

#define SIMD_SQRLENGTH( a, r )							\
					{SIMD_DOT_PRODUCT((a),(a),(r));}

#define SIMD_LENGTH( a, r )								\
					{SIMDVector sqrlen;					\
					 SIMD_SQRLENGTH((a),sqrlen);		\
					(r) = _mm_sqrt_ps(sqrlen);}

#define SIMD_NORMALIZE( a, r )											\
					{SIMDVector len;										\
					 SIMD_LENGTH((a),len);									\
					 len = _mm_shuffle_ps(len,len,_MM_SHUFFLE(0,0,0,0));	\
					 (r) = _mm_div_ps((a),len);	}

#define SIMD_NORMALIZE_FAST( a, r )											\
					{SIMDVector sqrlen;										\
					 SIMD_SQRLENGTH((a),sqrlen);							\
					 sqrlen = _mm_shuffle_ps(sqrlen,sqrlen,_MM_SHUFFLE(0,0,0,0));	\
					 (r) = _mm_mul_ps((a), _mm_rsqrt_ps(sqrlen) );	}

// // DO NOT USE IT : this is not precise at all ! >.<
//#define SIMD_NORMALIZE( a, r )										\
//					{SIMDVector l;									\
//					 SIMD_SQRLENGTH((a),l);							\
//					 l = _mm_rsqrt_ps(l);							\
//					 l = _mm_shuffle_ps(l,l,_MM_SHUFFLE(0,0,0,0));	\
//					 (r) = _mm_mul_ps((a),l);	}

#define SIMD_LERP( a, b, t, r )								\
					{ SIMDVector c = _mm_sub_ps((b),(a));	\
				      c = _mm_mul_ps(c,(t));				\
					 (r)= _mm_add_ps((a),c);}

//***************************************************************************
#define SIMD_MATRIX34_LOAD_A_12FLOAT(	float1_11, float1_12, float1_13, float1_14,				\
										float1_21, float1_22, float1_23, float1_24,				\
										float1_31, float1_32, float1_33, float1_34, r )			\
					{SIMD_LOAD_A_4FLOAT((float1_11), (float1_12), (float1_13), (float1_14), (r)->v[0] );	\
					 SIMD_LOAD_A_4FLOAT((float1_21), (float1_22), (float1_23), (float1_24), (r)->v[1] );	\
					 SIMD_LOAD_A_4FLOAT((float1_31), (float1_32), (float1_33), (float1_34), (r)->v[2] );}

#define SIMD_MATRIX34_ADD( m, n, r )							\
					{(r)->v[0] = _mm_add_ps( (m)->v[0], (n)->v[0] );	\
					 (r)->v[1] = _mm_add_ps( (m)->v[1], (n)->v[1] );	\
					 (r)->v[2] = _mm_add_ps( (m)->v[2], (n)->v[2] );}

#define SIMD_MATRIX34_MULT_SCALAR( m, v, r )				\
					{(r)->v[0] = _mm_mul_ps( (m)->v[0], (v) );	\
					 (r)->v[1] = _mm_mul_ps( (m)->v[1], (v) );	\
					 (r)->v[2] = _mm_mul_ps( (m)->v[2], (v) );}
 
  
#define SIMD_MATRIX34_MULT_VEC4( m, v, r )							\
					{SIMDVector a,b,c;								\
					 SIMD_DOT_PRODUCT((v),(m)->v[0],a);				\
					 SIMD_DOT_PRODUCT((v),(m)->v[1],b);				\
					 SIMD_DOT_PRODUCT((v),(m)->v[2],c);				\
					 a = _mm_unpacklo_ps(a,b);						\
					 (r) = _mm_shuffle_ps(a,c,_MM_SHUFFLE(1,0,1,0));}

		
typedef	__m128		SIMDVector;




		/// Matrix typedef !!! TODO : create nifty Matrix44 class, and use the VMX instructrions on XBOX360 !!
		//typedef XMMATRIX Matrix44;
		ALIGN_16( struct Matrix44 )
		{
			union
			{ 
				struct
				{
					REAL _11, _12, _13, _14;
					REAL _21, _22, _23, _24;
					REAL _31, _32, _33, _34;
					REAL _41, _42, _43, _44;
				};
				REAL m[4][4];
				REAL ma[16];
			};

			// static functions ---------------------------------
			inline static Matrix44*	Allocate( UINT32 _ulCount );
			inline static void		Free	( Matrix44* _pData );
		};


		inline Matrix44*	AllocateMatrix44( UINT32 _ulCount );
		inline void			FreeMatrix44	( Matrix44* _pData );


		ALIGN_16( struct Matrix43 )
		{
			union
			{ 
				struct
				{
					REAL _11, _12, _13;
					REAL _21, _22, _23;
					REAL _31, _32, _33;
					REAL _41, _42, _43;
				};
				REAL m[4][3]; 
			};

			// static functions ---------------------------------
			inline static Matrix43*	Allocate( UINT32 _ulCount );
			inline static void		Free	( Matrix43* _pData );
		};

		ALIGN_16( struct SIMDMatrix34 )
		{
			union
			{ 
				SIMDVector	v[3];
				//REAL		m[3][4]; 
				//struct
				//{
				//	REAL _11, _12, _13, _14;
				//	REAL _21, _22, _23, _24;
				//	REAL _31, _32, _33, _34;
 			//	};
			};

			// static functions ---------------------------------
			//static Matrix43*	Allocate( UINT32 _ulCount );
			//static void			Free	( Matrix43* _pData );

			inline void Load ( const Matrix43& _rMat )
			{
				SIMD_LOAD_A_4FLOAT( _rMat.m[0][0], _rMat.m[1][0], _rMat.m[2][0], _rMat.m[3][0], v[0] );
				SIMD_LOAD_A_4FLOAT( _rMat.m[0][1], _rMat.m[1][1], _rMat.m[2][1], _rMat.m[3][1], v[1] );
				SIMD_LOAD_A_4FLOAT( _rMat.m[0][2], _rMat.m[1][2], _rMat.m[2][2], _rMat.m[3][2], v[2] );
			}

			//void Store( Matrix43& _rMat )
			//{
			//	StoreFourFloats( & v[0], & _rMat.m[0][0], & _rMat.m[1][0], & _rMat.m[2][0], & _rMat.m[3][0] );
			//	StoreFourFloats( & v[1], & _rMat.m[0][1], & _rMat.m[1][1], & _rMat.m[2][1], & _rMat.m[3][1] );
			//	StoreFourFloats( & v[2], & _rMat.m[0][2], & _rMat.m[1][2], & _rMat.m[2][2], & _rMat.m[3][2] );
			//}
		};

		struct Matrix33
		{
			union
			{ 
				struct
				{
					REAL _11, _12, _13;
					REAL _21, _22, _23;
					REAL _31, _32, _33;
 				};
				REAL m[3][3];
				REAL ma[9];
			};
		};
		 
		inline bool Mtx_Equal( const Matrix44& _rM1, const Matrix44& _rM2 );
		
		inline void Mtx_BuildNull( Matrix44& matrixResult ); 
		inline void Mtx_BuildIdentity( Matrix44& matrix );
		
		inline void Mtx_BuildTranslation( Matrix44& _rOutMatrix, const Vec3& vPos );
		inline void Mtx_BuildScale( Matrix44& _rMatrix, const Vec3& vScale );
		inline void Mtx_BuildRotationX( Matrix44& _rOutMatrix, REAL _fAngleRadian );
		inline void Mtx_BuildRotationY( Matrix44& _rOutMatrix, REAL _fAngleRadian );
		inline void Mtx_BuildRotationZ( Matrix44& _rOutMatrix, REAL _fAngleRadian );
		inline void Mtx_BuildFromSRT( Matrix44& _rMatrix, const Vec3& _vTranslation, const Quaternion & _quatRotation, const Vec3 & _vScale);
		inline void Mtx_BuildFromSRT( Matrix44& _rMatrix, const SRT& _rSRT );

		inline void Mtx_BuildRotation( Matrix44& _rOutMatrix, const Vec3& _rvAxis, REAL _fAngleRadian );
		inline void Mtx_BuildPerspectiveRH( Matrix44& _rOutResult, REAL _fFovAngleY, REAL _fAspectRatio, REAL _fNearZ, REAL _fFarZ );
		inline void Mtx_BuildOrtographicRH( Matrix44& _rOutResult, REAL _fViewWidth, REAL _fViewHeight, REAL _fNearZ, REAL _fFarZ );

		inline void Mtx_Concat( const Matrix44& matrixOne, const Matrix44& matrixTwo, Matrix44& matrixResult ); 
		inline void Mtx_Transpose( const Matrix44& _rMatrix, Matrix44& _rOutResult ); 
		inline void Mtx_Inverse( const Matrix44& _rMatrix, Matrix44& _rOutResult , Vec3& _rvOutDeterminant );
		inline void Mtx_Scale( Matrix44& _rOutMatrix, const Matrix44& _rBaseMatrix, REAL _fScale );
		inline void Mtx_MultVector( const Matrix44& _rMatrix, const Vec3& _rvVec, Vec3& _rvOutResult );

  
		inline void Mtx_Add( const Matrix44& matrixOne, const Matrix44& matrixTwo, Matrix44& matrixResult ); 


		inline void Mtx_BuildFromSRT( Matrix43& _rMatrix, const SRT& _rSRT );
		inline void Mtx_BuildNull( Matrix43& matrixResult );		
		inline void Mtx_Scale( Matrix43& _rOutMatrix, const Matrix43& _rBaseMatrix, REAL _fScale );
		inline void Mtx_Add( const Matrix43& matrixOne, const Matrix43& matrixTwo, Matrix43& matrixResult ); 
		inline void Mtx_MultVector( const Matrix43& matrix, const Vec3& vec, Vec3& vecResult );

		inline void Mtx_MultVector( const Matrix33& matrix, const Vec3& vec, Vec3& vecResult );

		inline bool Mtx_UnitTest();


		inline void MatrixMulScalar			( const SIMDMatrix34& _pMatrix, REAL _fScalar, SIMDMatrix34& _pOutMatrix );
 		inline void MatrixAdd				( const SIMDMatrix34& _pMatA, const SIMDMatrix34& _pMatB, SIMDMatrix34& _pOutMatrix );
 		inline void MatrixMulVec			( const SIMDMatrix34& _pMatrixA, const Vec3& _pSrcVecA, Vec3& _pOutVecA );
 		inline void MatrixMulVecNormalize	( const SIMDMatrix34& _pMatrixA, const Vec3& _pSrcVecA, Vec3& _pOutVecA );
 		inline void MatrixMulVecNormalizeFast( const SIMDMatrix34& _pMatrixA, const Vec3& _pSrcVecA, Vec3& _pOutVecA );
 		inline void Vec3Normalize			( const Vec3& _rVec, Vec3& _rOutVec );
 		inline void Vec3NormalizeFast		( const Vec3& _rVec, Vec3& _rOutVec );
 
		inline bool SIMDUnitTest();
	}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	Declarations :	QUATERNION - quite useful to work and encode an orientation
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting
{
	namespace Math
	{ 
		/** Quaternion class, quite useful to work and encode an orientation 		*/
		class Quaternion
		{
		public:
			//---------------------------------------------------------------------

			/// CTOR - default
			inline Quaternion( void );
			/// CTOR
			inline Quaternion( REAL _fRealPart, REAL _fImaginaryPartX, REAL _fImaginaryPartY, REAL _fImaginaryPartZ );
			/// CTOR
			inline Quaternion( const Quaternion & _rQuat );

			/// CTOR
			inline Quaternion( Radian rAngle , const Vec3& _rvAxis );
			/// CTOR
			inline Quaternion( REAL _angleX, REAL _angleY, REAL _angleZ );

			//---------------------------------------------------------------------

			/// Construct a Quaternion from the given Euler Angles
			inline void					FromAngles		( const Vec3 & _rvAngles );
			/// Construct a Quaternion from the given Euler Angles
			inline void					FromAngles		( REAL _angleX, REAL _angleY, REAL _angleZ );
			/// Construct a Quaternion representing the shift from a direction to another
			inline void					FromDirections	( const Vec3 & _rvOldDirection, const Vec3 & _rvNewDirection );
			/// Construct a Quaternion representing orientation encoded in the given rotatin matrix - WARNING : determinant must be 1, so that no scale factor is encoded !
			inline void					FromMatrix		( const Matrix33& _rMatrix );

			/// Store the orientation information in the given Matrix
			inline void					ToMatrix		( Matrix44 & _rOutMatrix ) const;
			/// Store the orientation information in the given Matrix
			inline void					ToMatrix		( Matrix43 & _rOutMatrix ) const;
			//---------------------------------------------------------------------
			/// Rotate a vector by the current quaterion rotation
			inline Vec3					Rotate			( const Vec3 & _rvOldDirection ) const;

			//---------------------------------------------------------------------
			/// Get the REAL Part		( same as	'GetW()' )
			inline REAL					GetRealPart		( void ) const;
			/// Get the Imaginary Part	( same as	'Vec3( GetX(), GetY(), GetZ() )' )
			inline Vec3					GetImaginaryPart( void ) const;

			/// Get the X	- NB : same as 'GetImaginaryPart()[0]' 
			inline REAL					GetX			( void ) const;
			/// Get the Y	- NB : same as 'GetImaginaryPart()[1]' 
			inline REAL					GetY			( void ) const;
			/// Get the Z	- NB : same as 'GetImaginaryPart()[2]' 
			inline REAL					GetZ			( void ) const;
			/// Get the W	- NB : same as 'GetRealPart()'
			inline REAL					GetW			( void ) const;

			//---------------------------------------------------------------------	
			//
			inline bool				operator==		( const Quaternion & _rOther );
			//
			inline bool				operator!=		( const Quaternion & _rOther );

			inline Quaternion		operator-		( void ) const;

			inline Quaternion		operator+		( const Quaternion & _rQuat ) const;
			inline Quaternion		operator-		( const Quaternion & _rQuat ) const;	

			inline Quaternion		operator/		( REAL _fScalar ) const;	
			inline Quaternion		operator*		( REAL _fScalar ) const;

	        inline Vec3				operator* 		( const Vec3& _rvVec) const;
			inline Quaternion		operator*		( const Quaternion & _rQuat )	const;		


			inline Quaternion&		operator+=		( const Quaternion & _rQuat );
			inline Quaternion&		operator-=		( const Quaternion & _rQuat );

			inline Quaternion&		operator/=		( REAL _fScalar );
			inline Quaternion&		operator*=		( REAL _fScalar )	;

			inline Quaternion &		operator*=		( const Quaternion& _rQuat );


			//---------------------------------------------------------------------
			/// Inverse the Quaterion 
			inline Quaternion		GetInverse		( void ) const;
			/// Magnitude
			inline REAL				GetLength		( void ) const;

			/// Normalize the Quaternion 
			inline Quaternion &		Normalize		( void ); 

			/// Dot product 
			static inline REAL		DotP			( const Quaternion& _rQuatA, const Quaternion& _rQuatB ); 

 
			//---------------------------------------------------------------------

			/// Inverse quaternion, also known as conjugate
			inline static void		Inverse			( const Quaternion& _rBase, Quaternion& _rInverse );

			/// Spherical Linear Interpolation			
			inline static void		SLERP			( REAL _fTime , const Quaternion & _rStart, const Quaternion & _rEnd, Quaternion & _rOutResult ) ;

			/// Normalized Linear Interpolation			
			inline static void		NLERP			( REAL _fTime , const Quaternion & _rStart, const Quaternion & _rEnd, Quaternion & _rOutResult ) ;

		private:
			union 
			{
				struct
				{
					REAL	mfRealPart;			///< Real part
					Vec3	mvImaginaryPart;	///< Imaginary part
				};
				REAL ma[4];
			};
								
#ifndef __FINALROM		
		public:
			/// Print in the console	
			inline void 			Print( void ) const;
			/// Unit test the methods of this class
			inline static bool 		UnitTest( void );
#endif	
		};


	}// namespace Math
}// namespace Melting 


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	Declarations :	SRT - Scale Rotate Translate 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting 
{
	namespace Math
	{
		/** 
		 	SRT = Scale Rotate Translate 
		 */
		class SRT
		{ 
		public:
			
			/// Default CTOR : /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !
			inline SRT( void );
			
			/// CTOR - NDAL : I hate Euler angles, but since we have no Quaternions so far, we will stick to these
			inline SRT( const Vec3 & _rScale, const Vec3 & _rEulerAngles, const Vec3 & _rTranslate );	
			/// CTOR 
			inline SRT( const Vec3 & _rScale, const Quaternion & _rquatRotation, const Vec3 & _rTranslate );
			/// Copy CTOR
			inline SRT( const SRT & _rSRT ); 

			//
			inline bool operator==( const SRT & _rOther );
			//
			inline bool operator!=( const SRT & _rOther );

			
			
			/// Make the given SRT an Identity one (~Matrix Identity)
			static inline	void	BuildIdentity( SRT& _rOutIdentity );

			/// Make the given SRT an Null one (~every param is 0)
			static inline	void	BuildNull	( SRT& _rOutNull );

			/// Build the global SRT from the given _rParent and the given _rLocal, and stores the result in _rOutResult
			static inline	void	Transform	( const SRT& _rLocal, const SRT& _rParent, SRT& _rOutResult );

			/// Multiply the Base SRT by factor and store the result in _rOutResult
			static inline	void	Multiply		( REAL _fWeight, const SRT& _rBase, SRT& _rOutResult );
			/// Multiply the Base SRT by factor and store the result in _rOutResult
			static inline	void	MultiplyInPlace	( REAL _fWeight, SRT& _rOutResult );
			/// Sum up the given SRT and store the result in _rOutResult
			static inline	void	Add				( const SRT& _rA, const SRT& _rB, SRT& _rOutResult );
			/// Sum up the given SRT and store the result in _rOutResult
			static inline	void	AddInPlace		( const SRT& _rA, SRT& _rOutResult );

			/// Interpolate between _rStart abd _rEnd, at the given time (between 0.0f and 1.0f), and store the result in _rOutResult
			static inline	void	Interpolate	( REAL _fTime, const SRT& _rStart, const SRT& _rEnd, SRT& _rOutResult );
			
			/// Inverse 
			static inline	void	Inverse		( const SRT& _rBase, SRT& _rOutResult );

			 
			
			Vec3		mvScale; 			///< Scale
			Quaternion	mquatRotation;		///< Rotate - Quaternion !
			Vec3		mvTranslate; 		///< Translate
				
#ifndef __FINALROM		
		public:
			/// Print in the console	
			inline void 				Print( void ) const;
			/// Unit test the methods of this class
			inline static bool 			UnitTest( void );
#endif	
		}; 	 
	} // namespace Math
		
} // namespace Melting
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	Declarations :	AABBOX - Axis-Aligned Bounding Box, represented by its Min and Max vertex (i.e. Min = 'LEFT-BOTTOM-FRONT', Max = 'RIGHT-TOP-REAR')	
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting
{
	namespace Math
	{
		/**
			Axis-Aligned Bounding Box, represented by its Min and Max vertex (i.e. Min = 'LEFT-BOTTOM-FRONT', Max = 'RIGHT-TOP-REAR')	
		*/
		class AABBox
		{	
		
		public:
			
			
			/// Default CTOR : /!\ Leave it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !
			inline AABBox( void );
			
			/// CTOR
			inline AABBox( const Vec3 & _rvMin, const Vec3 & _rvMax );
			
			/// Copy-CTOR
			inline AABBox( const AABBox & _rBox );
			
						
			/// Setup the AABBox with the given 'min' position (bottom left front) of the box, and its volume (with absolute components required!)
			inline void 			SetupFromMinPositionAndVolume		( const Vec3 & _rvMinPosition, const Vec3 & _rvVolume );
			
			/// Setup the AABBox with the given 'center' position of the box	, and its volume / 2 (with absolute components required!)
			inline void 			SetupFromCenterPositionAndHalfVolume( const Vec3 & _rvCenterPosition, const Vec3 & _rvHalfVolume );
			

			/// Get the Min...
			inline const Vec3&		GetMin( void ) const;
			/// Get the Max !!!
			inline const Vec3&		GetMax( void ) const;

			/// Set the Min...
			inline void				SetMin( const Vec3& ) ;
			/// Set the Max !!!
			inline void				SetMax( const Vec3& ) ;

			/// Get the Center of the AABBox
			inline void 			GetCenter( Vec3& _rvResult ) const;
			
			/// Get Half the height (Y) of this AABBox
			inline REAL				GetHalfHeight( void ) const;
					
			/// Get the 'Positive vertex' according to the given normal
			inline void				GetPositiveVertex	( const Vec3& _rvNormal, Vec3& _rvResult ) const;
			/// Get the 'Negative vertex' according to the given normal
			inline void				GetNegativeVertex	( const Vec3& _rvNormal, Vec3& _rvResult ) const;			
			/// Snap the given Position to the closest plane of the AABBox
			inline void 			SnapToClosestPlane	( bool _abSkippedDimensions[3], Vec3 & _rvOutSnaped ) const;
			
			/// Get the center of the model, and store it inside _rvOutCenter
			inline void				ComputeCenter( Vec3&  _rvOutCenter ) const;
			
			/// Is the given point inside the AABBox
			inline bool				IsInside( const Vec3& _rvPoint ) const;
				
				
			/// operator : ==
			inline bool 			operator == ( const AABBox & _rBox ) const;			
			/// operator : !=
			inline bool 			operator != ( const AABBox & _rBox ) const;
		
		
			
			/// Compute the union of two AABBox, and store the result in the given reference parameter 'rOutBoxResult'
			inline static void 		Union( const AABBox & _rBoxA, const AABBox & _rBoxB, AABBox & _rOutBoxResult );

			/// Compute the union of a given AABBox and a point, and store the result in the given reference parameter 'rOutBoxResult'
			inline static void 		Union(const AABBox & _rBoxA, const Vec3 & _vPoint, AABBox & _rOutBoxResult);
			

		private:
				
			Vec3		mvMin;			///< Point 'LEFT-BOTTOM-FRONT'
			Vec3		mvMax;			///< Point 'RIGHT-TOP-REAR'	
		
#ifndef __FINALROM
		public:			
			/// Print infos to the console
			inline void 				Print( void ) const;
			/// Unit test the methods of this class
			inline static bool 			UnitTest( void );
#endif	
		}; 
	}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	Declarations :	BOUNDING SPHERE
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting
{
	namespace Math
	{ 			
		/** Simple Bounding Sphere 	*/
		class BoundingSphere
		{
		public:
		
			/// Default CTOR : /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !
			inline BoundingSphere( void );
			
			/// CTOR
			inline BoundingSphere( const Vec3& _rvCenter, REAL _fRadius ) ;
			
			/// Copy-CTOR
			inline BoundingSphere( const BoundingSphere& _rSphere ) ;
 
			/// Setup this BoundingSphere as the tightest sphere containing the given box
			inline void 			SetupFromAABBoxInside( const AABBox & rBox );
			
			/// Setup this BoundingSphere as the biggest sphere contained in the given box
			inline void 			SetupFromAABBoxOutside( const AABBox & rBox );
			
			/// Get the BoundingSphere center
			inline	const Vec3& 	GetCenter( void ) const;
			/// Set the BoundingSphere center
			inline	void 			SetCenter( const Vec3& _rvCenter );
			/// Get the BoundingSphere radius
			inline	REAL 			GetRadius( void ) const;
			/// Set the BoundingSphere radius
			inline	void  			SetRadius( REAL _fRadius ) ; 


			/// Is the given point inside the AABBox
			inline bool				IsInside( const Vec3& _rvPoint ) const;
			
			/// operator : ==	
			inline bool 			operator == ( const BoundingSphere& _rBoundingSphere ) const;	
			/// operator : !=	
			inline bool 			operator != ( const BoundingSphere& _rBoundingSphere ) const;
			
		
			/// Compute the union of rSphereA and rSphereB, and store the result in rOutUnionSphere -- PS : rSphereA or rSphereB can be used as rOutUnionSphere 
			inline static void 		Union( const BoundingSphere& _rSphereA , const BoundingSphere& _rSphereB, BoundingSphere& _rOutUnionSphere );
			
	
		private:
		
			Vec3 		mvCenter;			///< center of the sphere
			REAL		mfRadius;			///< Radius R 
		
		
#ifndef __FINALROM
		public:			
			/// Print infos to the console
			inline void 				Print( void ) const;
			/// Unit test the methods of this class
			inline static bool 			UnitTest( void );
#endif	
		};   
	}
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	Declarations :	CONE 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting
{
	namespace Math
	{ 			 
		/**
			A Cone is a finite space, represented by 
				- a vertex (its starting point)
				- a direction (unit vector)
				- a length along the direction 
				- an opening angle 
		*/
		class Cone
		{ 
		public: 
			
			/// Default CTOR : /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !		
			inline Cone( void )  ;
			
			/// CTOR
			inline Cone( const Vec3 & _rvVertex, const Vec3 & _rvDirection, REAL fLength, REAL _fAngle );
			
			/// Copy CTOR
			inline Cone( const Cone& _rCone );
			
			/// Setup the cone from a starting point (the vertex), and endind point, and an angle
			inline void 			SetupFromTwoPointsAndAngle( const Vec3 & _rvStart, const Vec3 & _rvEnd, REAL _fAngle );
			
			
			inline void				SetVertex		( const Vec3 & _rvVertex );
			inline const Vec3& 		GetVertex		( void ) const ;
			
			inline void				SetDirection	( const Vec3 & _rvDirection ) ;
			inline const Vec3& 		GetDirection	( void ) const ;
			
			inline void 			SetLength		( REAL fLength );
			inline REAL				GetLength		( void ) const;			
			
			inline void 			SetAngle		( REAL _fAngle );
			inline REAL				GetAngle		( void ) const;
			inline REAL				GetHalfAngleCos	( void ) const;
			inline REAL				GetHalfAngleSin	( void ) const;
			
			
			/// operator : == 
			inline bool				operator == ( const Cone & _rCone ) const;
			/// operator : !=  
			inline bool				operator != ( const Cone & _rCone ) const;
		
		
			/// Is the Sphere inside the Cone ??
			inline static bool		IsInside		( const BoundingSphere & _rSphere, const Cone& _rCone );
			
			/// Is the Point inside the Cone ??
			inline static bool		IsInside		( const Vec3& _rPoint, const Cone& _rCone );
			
		private:
			/// Compute the Sin & Cos of the (Angle/2) 
			inline void				ComputeSinCos	( void ); 
			
		
			Vec3		mvVertex; 			///< pic 
			Vec3		mvDirection;		///< unit-length direction
			REAL		mfLength;			///< length or height
			REAL		mfAngle;			///< angle
			
			REAL		mfHalfAngleCos;		///< cos( angle /2 ) - computed once and then stored
			REAL		mfHalfAngleSin;		///< sin( angle /2 ) - computed once and then stored
			
#ifndef __FINALROM		
		public:
			/// Print in the console	
			inline void 			Print( void ) const;
			/// Unit test the methods of this class
			inline static bool 		UnitTest( void );
#endif	
		}; 
	}
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	Declarations :	PLANE 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting
{
	namespace Math
	{
		/** 
		 	A plane is infinite, and simply defined by its normal, and a factor d telling "where this plane lies along its normal" 
		 */
		class Plane
		{ 
		public:
		
			/// Default CTOR : /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !
			inline 	Plane ( void );
			
			/// CTOR
			inline 	Plane ( const Vec3 & _rvNormal , REAL _fDistanceToOrigin ); 
			
			/// Copy - CTOR
			inline 	Plane ( const Plane & _rPlane );
		
		

			/// Create a Plane defined by a triangle
			inline 	void 			SetupFromTriangle		( const Vec3 & _rvVertex1, const Vec3 & _rvVertex2, const Vec3 & _rvVertex3 );
			/// Create a Plane defined by a normal and a point on this Plane
			inline 	void 			SetupFromNormalAndPoint	( const Vec3 & _rvNormal, const Vec3 & _rvVertex );
	
			/// Get the Plane normal -- NO SETTER ! use one of the method above to setup the Plane
			inline 	const Vec3& 	GetNormal 	( void ) const;
			/// Get the Plane distance to the world origin (along the direction) -- NO SETTER ! use one of the method above to setup the Plane
			inline 	REAL 			GetDistance	( void ) const;
			
			
			/// Compute the distance from the point to this Plane 
			inline 	REAL			DistanceFrom( const Vec3 & _rvVertex ) const;

			/// Project the given point 'rvVertex 'on this Plane, and store the _result in 'rvOutProjectedVertex';
			inline 	void 			Project		( const Vec3 & _rvVertex , Vec3 & _rvOutProjectedVertex );
		
		
			///  Tells wether the given point is above this Plane -- /!\ Assume the normal points toward the "ABOVE"  or "FRONT" 
			inline 	bool 			IsAbove		( const Vec3 & _rvVertex ) const; 


			/// operator : == 
			inline 	bool 			operator == ( const Plane & _rPlane ) const;
			
			/// operator : != 
			inline 	bool 			operator != ( const Plane & _rPlane ) const;
			
		private:			
			
			Vec3	mvNormal;		///< Normal
			REAL	mfDist; 		///< Distance from the Plane to the center of the world
		
		
#ifndef __FINALROM	
		public:	
			/// Print in the console	
			inline void 				Print( void ) const;
			/// Unit test the methods of this class
			inline static bool 			UnitTest( void );
#endif	 

		}; 	 
	} // namespace Math
		
} // namespace Melting
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//	Declarations :	RAY 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting
{
	namespace Math
	{
		/** 
		 	A Ray is infinite, and simply defined by its direction (so it's oriented), and a point (let's say it's its origin)
		 */
		class Ray
		{ 
		public:
			
			/// Default CTOR : /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !
			inline Ray( void );
			
			/// CTOR
			inline Ray( const Vec3 & _rvOrigin, const Vec3 & _rvDirection );
			
			/// Copy CTOR
			inline Ray( const Ray & _rRay );
		
			/// Create a Ray defined by a segment composed of 2 points : Start & End
			inline void 				SetupFromSegment	( const Vec3 & _rvVertexStart, const Vec3 & _rvVertexEnd );
			
			/// Set the origin of the Ray
			inline void					SetOrigin			( const Vec3 & _rvOrigin ) ;
			/// Get the origin of the Ray
			inline const Vec3 & 		GetOrigin			( void ) const;
			/// Set the direction of the Ray
			inline void			 		SetDirection		( const Vec3 & _rvDirection );
			/// Get the direction of the Ray
			inline const Vec3 & 		GetDirection		( void ) const;
			/// Get the 1/direction of the Ray
			inline const Vec3 & 		GetInverseDirection	( void ) const ;
			
			/// Compute the 3D position of a point on the Ray, located at the given distance from the Ray origin -- the _result is stored in the _reference parameter 'rvOutPoint'
			inline void 				ComputePoint		( REAL _fDistanceFromOrigin, Vec3 & _rvOutPoint ) const;
			/// Compute the distance from the Ray origin to the given Plane
			inline REAL					ComputeDistance		( const Plane & _rPlane ) const;
			
			/// Compute the intersection of the given Plane with this Ray , and store the _result in 'rvOutIntersect', and _returns wether there was an intersection or not
			inline bool 				ComputeIntersection	( const Plane & _rPlane, Vec3 & _rvOutIntersect ) const;
			/// Returns wether this Ray intersects the given Plane or not
			inline bool 				DoesIntersect		( const Plane & _rPlane ) const;
		 
			/// Returns wether this Ray intersects the given AABBox or not
			inline bool 				DoesIntersect		( const AABBox & _rAABBox ) const;
			
			/// Intersect _ray R(t) = p + t*d against AABB a. When intersecting,
			/// _return intersection distance tmin and point q of intersection
			inline bool 				ComputeIntersection ( const AABBox & _rAABBox, Vec3& _rvOutIntersection, REAL& _rOutIntersectionDistance ) const;
			
			/// same as ComputeIntersection, except that we can specify the ScaleRotateTranslate of the box !
			inline bool 				ComputeIntersection ( const AABBox & _rAABBox, const SRT& _rBoxSRT, Vec3& _rvOutIntersection, REAL& _rOutIntersectionDistance ) const;
 		
			/// Returns wether this Ray intersects the given BoundingSphere or not
			inline bool 				DoesIntersect		( const BoundingSphere & _rBSphere ) const;
	
	
			/// operator : ==
			inline bool 				operator == ( const Ray & _rRay ) const;
			/// operator : !=
			inline bool 				operator != ( const Ray & _rRay ) const;
			
	
		private : 			
			
			Vec3		mvOrigin; 			///< Origin of the Ray
			Vec3		mvDirection;		///< Direction in which the Ray in going	
			Vec3		mvInverseDirection;	///< 1/Direction
			Vec3UINT8	mvSignedNeg;		///< Sign of the Ray direction ( 1 == negative, 0 == positive )
					
#ifndef __FINALROM		
		public:
			/// Print in the console	
			inline void 				Print( void ) const;
			/// Unit test the methods of this class
			inline static bool 			UnitTest( void );
#endif	
		}; 	 
	} // namespace Math
		
} // namespace Melting




		


//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//	Implementations :	UTILS 
//'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 
namespace Melting
{
	namespace Math
	{
		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		bool Equal( REAL _realA, REAL _realB, REAL _epsilon )
		{
			REAL fAbs = Math::Abs( _realA - _realB );
			bool bEqual = fAbs < _epsilon;
			
			return Math::Abs( _realA - _realB ) < _epsilon;
		}
		// ------------------------------------------------------------------------------ 
		// ------------------------------------------------------------------------------
		REAL Inverse( REAL _real )
		{
			return REAL( 1.0f ) /_real;
		}
		
		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		void Swap( REAL& _rVal1, REAL& _rVal2 )
		{
			REAL Tmp( _rVal1 );			
			_rVal1 = _rVal2;
			_rVal2 = Tmp;
		} 
		
		// ------------------------------------------------------------------------------ 
		// ------------------------------------------------------------------------------
		void Clamp( REAL & _rfOutClamped, REAL _fMin, REAL _fMax )
		{
			MB_ASSERT( _fMin < _fMax, "Min and Max inverted !" );

			if( _rfOutClamped < _fMin )
			{
				_rfOutClamped = _fMin;
			}
			else
			if( _fMax < _rfOutClamped )
			{
				_rfOutClamped = _fMax;
			}
		}

		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		Radian DegreeToRadian( Degree _fAngle )
		{
			return _fAngle * PI / 180.0f;
		}

		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		Degree RadianToDegree( Radian _fAngle )
		{
			return _fAngle * 180.0f / PI ;
		}

		// ------------------------------------------------------------------------------ 
		// ------------------------------------------------------------------------------
		REAL Sin( Radian _fAngle )
		{
			return sin( _fAngle ); 
		}	
		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		Radian ASin( REAL _fValue )
		{
			return asin( _fValue ); 
		}

		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		REAL Cos( Radian _fAngle )
		{ 
			return cos( _fAngle );
		}
 
		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		Radian ACos( REAL _fValue )
		{
			return acos( _fValue ) ;
		} 

		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		REAL Tan( Radian _fAngle )
		{ 
			return tan( _fAngle );
		}

		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		Radian ATan( REAL _fValue )
		{
			return atan( _fValue ) ;
		}
		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		Radian ATan2( REAL _fValueY, REAL _fValueX )
		{
			return atan2( _fValueY, _fValueX ) ;
		}

		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		REAL Square( REAL _real )
		{
			return _real*_real;
		}
		
		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		REAL SquareRoot( REAL _real )
		{
			return sqrt( _real );
		}
		
		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		REAL InverseSquareRoot( REAL _real )
		{
			return Inverse( SquareRoot ( _real ) );
		}

		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		REAL Abs( REAL _real )
		{
			return fabs(_real);
		}
		
		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		UINT32 NextMultiple( UINT32 _ulValue, UINT32 _ulMultiple )
		{
			return ((_ulValue + _ulMultiple - 1) / _ulMultiple) * _ulMultiple;
		}
	}
}


//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//	Implementations :	VECTOR

// TODO .......................

//'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 

//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//	Implementations :	MATRIX 
//'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 
namespace Melting
{
	namespace Math
	{		
		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		void Mtx_Scale( Matrix43& _rOutMatrix, const Matrix43& _rBaseMatrix, REAL _fScale )
		{	  
			_rOutMatrix.m[0][0] = _rBaseMatrix.m[0][0] * _fScale;
			_rOutMatrix.m[0][1] = _rBaseMatrix.m[0][1] * _fScale;
			_rOutMatrix.m[0][2] = _rBaseMatrix.m[0][2] * _fScale;

			_rOutMatrix.m[1][0] = _rBaseMatrix.m[1][0] * _fScale;
			_rOutMatrix.m[1][1] = _rBaseMatrix.m[1][1] * _fScale;
			_rOutMatrix.m[1][2] = _rBaseMatrix.m[1][2] * _fScale;

			_rOutMatrix.m[2][0] = _rBaseMatrix.m[2][0] * _fScale;
			_rOutMatrix.m[2][1] = _rBaseMatrix.m[2][1] * _fScale;
			_rOutMatrix.m[2][2] = _rBaseMatrix.m[2][2] * _fScale;

			_rOutMatrix.m[3][0] = _rBaseMatrix.m[3][0] * _fScale;
			_rOutMatrix.m[3][1] = _rBaseMatrix.m[3][1] * _fScale;
			_rOutMatrix.m[3][2] = _rBaseMatrix.m[3][2] * _fScale;
 		}

		void Mtx_MultVector( const Matrix43& _rMatrix, const Vec3& vec, Vec3& vecResult )
		{
 			vecResult[ 0 ] =	vec[ 0 ] * _rMatrix.m[ 0 ][ 0 ] 
							+	vec[ 1 ] * _rMatrix.m[ 1 ][ 0 ] 
							+	vec[ 2 ] * _rMatrix.m[ 2 ][ 0 ] 
							+	/* 1* */ _rMatrix.m[ 3 ][ 0 ];

			vecResult[ 1 ] =	vec[ 0 ] * _rMatrix.m[ 0 ][ 1 ] 
							+	vec[ 1 ] * _rMatrix.m[ 1 ][ 1 ] 
							+	vec[ 2 ] * _rMatrix.m[ 2 ][ 1 ] 
							+	/* 1* */  _rMatrix.m[ 3 ][ 1 ];

			vecResult[ 2 ] =	vec[ 0 ] * _rMatrix.m[ 0 ][ 2 ] 
							+	vec[ 1 ] * _rMatrix.m[ 1 ][ 2 ] 
							+	vec[ 2 ] * _rMatrix.m[ 2 ][ 2 ] 
							+	/* 1* */  _rMatrix.m[ 3 ][ 2 ];
    	}

		void Mtx_Add( const Matrix43& matrixOne, const Matrix43& matrixTwo, Matrix43& matrixResult )
		{
  			matrixResult.m[0][0] = matrixOne.m[0][0] + matrixTwo.m[0][0];
			matrixResult.m[0][1] = matrixOne.m[0][1] + matrixTwo.m[0][1];
			matrixResult.m[0][2] = matrixOne.m[0][2] + matrixTwo.m[0][2];
 
			matrixResult.m[1][0] = matrixOne.m[1][0] + matrixTwo.m[1][0];
			matrixResult.m[1][1] = matrixOne.m[1][1] + matrixTwo.m[1][1];
			matrixResult.m[1][2] = matrixOne.m[1][2] + matrixTwo.m[1][2];
 
			matrixResult.m[2][0] = matrixOne.m[2][0] + matrixTwo.m[2][0];
			matrixResult.m[2][1] = matrixOne.m[2][1] + matrixTwo.m[2][1];
			matrixResult.m[2][2] = matrixOne.m[2][2] + matrixTwo.m[2][2];
 
			matrixResult.m[3][0] = matrixOne.m[3][0] + matrixTwo.m[3][0];
			matrixResult.m[3][1] = matrixOne.m[3][1] + matrixTwo.m[3][1];
			matrixResult.m[3][2] = matrixOne.m[3][2] + matrixTwo.m[3][2];
  		}
 
 
		void MatrixMulScalar( const SIMDMatrix34& _pMatrix, REAL _fScalar, SIMDMatrix34& _pOutMatrix )
		{
			SIMDVector v;
			SIMD_LOAD_A_FLOAT_AS_XXXX( _fScalar, v );
			
			SIMD_MATRIX34_MULT_SCALAR(&_pMatrix,v,&_pOutMatrix);
		}

		void MatrixAdd( const SIMDMatrix34& _pMatA, const SIMDMatrix34& _pMatB, SIMDMatrix34& _pOutMatrix )
		{
			SIMD_MATRIX34_ADD(&_pMatA,&_pMatB,&_pOutMatrix);
		}

		void MatrixMulVec(	const SIMDMatrix34& _pMatrixA, const Vec3& _pSrcVecA, Vec3& _pOutVecA )
		{ 
			SIMDVector v;
			SIMD_LOAD_A_FLOAT3_AS_1ZYX( _pSrcVecA, v );
									
			SIMDVector r;
			SIMD_MATRIX34_MULT_VEC4( &_pMatrixA, v, r ); 		
			
			SIMD_STORE_A_FLOAT3( _pOutVecA, r ); 
		} 

		
		void MatrixMulVecNormalize(	const SIMDMatrix34& _pMatrixA, const Vec3& _pSrcVecA, Vec3& _pOutVecA )
		{ 
			SIMDVector v;
			SIMD_LOAD_A_FLOAT3_AS_1ZYX( _pSrcVecA, v );
									
			SIMDVector r;
			SIMD_MATRIX34_MULT_VEC4( &_pMatrixA, v, r ); 		
			SIMD_NORMALIZE( r, v );
			
			SIMD_STORE_A_FLOAT3( _pOutVecA, v ); 
		} 

		void MatrixMulVecNormalizeFast(	const SIMDMatrix34& _pMatrixA, const Vec3& _pSrcVecA, Vec3& _pOutVecA )
		{ 
			SIMDVector v;
			SIMD_LOAD_A_FLOAT3_AS_1ZYX( _pSrcVecA, v );
									
			SIMDVector r;
			SIMD_MATRIX34_MULT_VEC4( &_pMatrixA, v, r ); 		
			SIMD_NORMALIZE_FAST( r, v );
			
			SIMD_STORE_A_FLOAT3( _pOutVecA, v ); 
		} 

		void Vec3Normalize( const Vec3& _rVec, Vec3& _rOutVec )
		{ 
			SIMDVector a;
			SIMD_LOAD_A_FLOAT3_AS_0ZYX( _rVec, a ); 
			  
			SIMDVector result;
			SIMD_NORMALIZE( a, result );

			SIMD_STORE_A_FLOAT3( _rOutVec, result ); 
		}

		void Vec3NormalizeFast( const Vec3& _rVec, Vec3& _rOutVec )
		{ 
			SIMDVector a;
			SIMD_LOAD_A_FLOAT3_AS_0ZYX( _rVec, a ); 
			  
			SIMDVector result;
			SIMD_NORMALIZE_FAST( a, result );

			SIMD_STORE_A_FLOAT3( _rOutVec, result ); 
		}
				
		inline UINT8* AlignedMalloc( UINT32 _iBlockSize, UINT32 _iAlignment )
		{
			UINT8* pMemory = (UINT8*) malloc( _iBlockSize + _iAlignment );

			UINT8* pAlignedMemory = (UINT8*)((UINT32)( pMemory + _iAlignment ) & ( ~ (_iAlignment - 1 ) ));

			// safe, because adding _iAlignment instead of (_iAlignment-1) guarantees it is past the actual start of the maloc'ed block.
			*(pAlignedMemory - 1) = pAlignedMemory - pMemory;

			return pAlignedMemory;
		}

		inline void AlignedFree( UINT8* _pAlignedMemory )
		{
			UINT8 iOffset = *(_pAlignedMemory - 1);
			free( (void*)(_pAlignedMemory - iOffset) );
		}
		 
		Matrix44* AllocateMatrix44( UINT32 _ulCount )
		{
			return (Matrix44*)	AlignedMalloc( _ulCount * sizeof(Matrix44), 16 );
		}

		void FreeMatrix44( Matrix44* _pData )
		{
			AlignedFree( (UINT8*)_pData );
		}  
  
		Matrix43* Matrix43::Allocate( UINT32 _ulCount )
		{
			//return (Matrix43*)	_aligned_malloc( _ulCount * sizeof(Matrix43), 16 );
			return (Matrix43*)	AlignedMalloc( _ulCount * sizeof(Matrix43), 16 );
		}

		void Matrix43::Free( Matrix43* _pData )
		{
			//_aligned_free( (void*)_pData );
			AlignedFree( (UINT8*)_pData );
		}  

		bool Mtx_Equal( const Matrix44& _rM1, const Matrix44& _rM2 )
		{
//			return * _rM1.ma == * _rM2.ma;

			return		Math::Equal(_rM1.ma[0], _rM2.ma[0] )
					&&	Math::Equal(_rM1.ma[1], _rM2.ma[1] )
					&&	Math::Equal(_rM1.ma[2], _rM2.ma[2] )
					&&	Math::Equal(_rM1.ma[3], _rM2.ma[3] )
					&&	Math::Equal(_rM1.ma[4], _rM2.ma[4] )
					&&	Math::Equal(_rM1.ma[5], _rM2.ma[5] )
					&&	Math::Equal(_rM1.ma[6], _rM2.ma[6] )
					&&	Math::Equal(_rM1.ma[7], _rM2.ma[7] )
					&&	Math::Equal(_rM1.ma[8], _rM2.ma[8] )
					&&	Math::Equal(_rM1.ma[9], _rM2.ma[9] )
					&&	Math::Equal(_rM1.ma[10], _rM2.ma[10] )
					&&	Math::Equal(_rM1.ma[11], _rM2.ma[11] )
					&&	Math::Equal(_rM1.ma[12], _rM2.ma[12] )
					&&	Math::Equal(_rM1.ma[13], _rM2.ma[13] )
					&&	Math::Equal(_rM1.ma[14], _rM2.ma[14] )
					&&	Math::Equal(_rM1.ma[15], _rM2.ma[15] )
					;
		}

			
		void Mtx_BuildNull( Matrix44& _rOutMatrix )
		{
 			_rOutMatrix.m[0][0] = 0.0f;
			_rOutMatrix.m[0][1] = 0.0f;
			_rOutMatrix.m[0][2] = 0.0f;
			_rOutMatrix.m[0][3] = 0.0f;

			_rOutMatrix.m[1][0] = 0.0f;
			_rOutMatrix.m[1][1] = 0.0f;
			_rOutMatrix.m[1][2] = 0.0f;
			_rOutMatrix.m[1][3] = 0.0f;

			_rOutMatrix.m[2][0] = 0.0f;
			_rOutMatrix.m[2][1] = 0.0f;
			_rOutMatrix.m[2][2] = 0.0f;
			_rOutMatrix.m[2][3] = 0.0f;

			_rOutMatrix.m[3][0] = 0.0f;
			_rOutMatrix.m[3][1] = 0.0f;
			_rOutMatrix.m[3][2] = 0.0f;
			_rOutMatrix.m[3][3] = 0.0f;
		}


		void Mtx_BuildIdentity( Matrix44& _rOutMatrix )
		{ 
			_rOutMatrix.m[0][0] = 1.0f;
			_rOutMatrix.m[0][1] = 0.0f;
			_rOutMatrix.m[0][2] = 0.0f;
			_rOutMatrix.m[0][3] = 0.0f;
			
			_rOutMatrix.m[1][0] = 0.0f;
			_rOutMatrix.m[1][1] = 1.0f;
			_rOutMatrix.m[1][2] = 0.0f;
			_rOutMatrix.m[1][3] = 0.0f;
			
			_rOutMatrix.m[2][0] = 0.0f;
			_rOutMatrix.m[2][1] = 0.0f;
			_rOutMatrix.m[2][2] = 1.0f;
			_rOutMatrix.m[2][3] = 0.0f;

			_rOutMatrix.m[3][0] = 0.0f;
			_rOutMatrix.m[3][1] = 0.0f;
			_rOutMatrix.m[3][2] = 0.0f;
			_rOutMatrix.m[3][3] = 1.0f;
		}

		void Mtx_BuildTranslation( Matrix44& _rOutM, const Vec3& _rvPos )
		{
			_rOutM.m[0][0] = 1.0f;
			_rOutM.m[0][1] = 0.0f;
			_rOutM.m[0][2] = 0.0f;
			_rOutM.m[0][3] = 0.0f;

			_rOutM.m[1][0] = 0.0f;
			_rOutM.m[1][1] = 1.0f;
			_rOutM.m[1][2] = 0.0f;
			_rOutM.m[1][3] = 0.0f;

			_rOutM.m[2][0] = 0.0f;
			_rOutM.m[2][1] = 0.0f;
			_rOutM.m[2][2] = 1.0f;
			_rOutM.m[2][3] = 0.0f;

			_rOutM.m[3][0] = _rvPos[0];
			_rOutM.m[3][1] = _rvPos[1];
			_rOutM.m[3][2] = _rvPos[2];
			_rOutM.m[3][3] = 1.0f;
		}

		void Mtx_BuildScale( Matrix44& _rOutM, const Vec3& _rvScale )
		{			
			_rOutM.m[0][0] = _rvScale[ 0 ];
			_rOutM.m[0][1] = 0.0f;
			_rOutM.m[0][2] = 0.0f;
			_rOutM.m[0][3] = 0.0f;

			_rOutM.m[1][0] = 0.0f;
			_rOutM.m[1][1] = _rvScale[ 1 ];
			_rOutM.m[1][2] = 0.0f;
			_rOutM.m[1][3] = 0.0f;

			_rOutM.m[2][0] = 0.0f;
			_rOutM.m[2][1] = 0.0f;
			_rOutM.m[2][2] = _rvScale[ 2 ];
			_rOutM.m[2][3] = 0.0f;

			_rOutM.m[3][0] = 0.0f;
			_rOutM.m[3][1] = 0.0f;
			_rOutM.m[3][2] = 0.0f;
			_rOutM.m[3][3] = 1.0f;
		} 

		void Mtx_Scale( Matrix44& _rOutMatrix, const Matrix44& _rBaseMatrix, REAL _fScale )
		{			 
 			_rOutMatrix.m[0][0] = _rBaseMatrix.m[0][0] * _fScale;
			_rOutMatrix.m[0][1] = _rBaseMatrix.m[0][1] * _fScale;
			_rOutMatrix.m[0][2] = _rBaseMatrix.m[0][2] * _fScale;
			_rOutMatrix.m[0][3] = _rBaseMatrix.m[0][3] * _fScale;

			_rOutMatrix.m[1][0] = _rBaseMatrix.m[1][0] * _fScale;
			_rOutMatrix.m[1][1] = _rBaseMatrix.m[1][1] * _fScale;
			_rOutMatrix.m[1][2] = _rBaseMatrix.m[1][2] * _fScale;
			_rOutMatrix.m[1][3] = _rBaseMatrix.m[1][3] * _fScale;

			_rOutMatrix.m[2][0] = _rBaseMatrix.m[2][0] * _fScale;
			_rOutMatrix.m[2][1] = _rBaseMatrix.m[2][1] * _fScale;
			_rOutMatrix.m[2][2] = _rBaseMatrix.m[2][2] * _fScale;
			_rOutMatrix.m[2][3] = _rBaseMatrix.m[2][3] * _fScale;

			_rOutMatrix.m[3][0] = _rBaseMatrix.m[3][0] * _fScale;
			_rOutMatrix.m[3][1] = _rBaseMatrix.m[3][1] * _fScale;
			_rOutMatrix.m[3][2] = _rBaseMatrix.m[3][2] * _fScale;
			_rOutMatrix.m[3][3] = _rBaseMatrix.m[3][3] * _fScale;
		}

		void Mtx_BuildRotationX( Matrix44& _rOutMatrix, REAL _fAngleRadian )
		{
			REAL fSinAngle = Math::Sin( _fAngleRadian );
			REAL fCosAngle = Math::Cos( _fAngleRadian );

			_rOutMatrix.m[0][0] = 1.0f;
			_rOutMatrix.m[0][1] = 0.0f;
			_rOutMatrix.m[0][2] = 0.0f;
			_rOutMatrix.m[0][3] = 0.0f;

			_rOutMatrix.m[1][0] = 0.0f;
			_rOutMatrix.m[1][1] = fCosAngle;
			_rOutMatrix.m[1][2] = fSinAngle;
			_rOutMatrix.m[1][3] = 0.0f;

			_rOutMatrix.m[2][0] = 0.0f;
			_rOutMatrix.m[2][1] = -fSinAngle;
			_rOutMatrix.m[2][2] = fCosAngle;
			_rOutMatrix.m[2][3] = 0.0f;

			_rOutMatrix.m[3][0] = 0.0f;
			_rOutMatrix.m[3][1] = 0.0f;
			_rOutMatrix.m[3][2] = 0.0f;
			_rOutMatrix.m[3][3] = 1.0f;
		}

		void Mtx_BuildRotationY( Matrix44& _rOutMatrix, REAL _fAngleRadian )
		{
			REAL fSinAngle = Math::Sin( _fAngleRadian );
			REAL fCosAngle = Math::Cos( _fAngleRadian );

			_rOutMatrix.m[0][0] = fCosAngle;
			_rOutMatrix.m[0][1] = 0.0f;
			_rOutMatrix.m[0][2] = -fSinAngle;
			_rOutMatrix.m[0][3] = 0.0f;

			_rOutMatrix.m[1][0] = 0.0f;
			_rOutMatrix.m[1][1] = 1.0f;
			_rOutMatrix.m[1][2] = 0.0f;
			_rOutMatrix.m[1][3] = 0.0f;

			_rOutMatrix.m[2][0] = fSinAngle;
			_rOutMatrix.m[2][1] = 0.0f;
			_rOutMatrix.m[2][2] = fCosAngle;
			_rOutMatrix.m[2][3] = 0.0f;

			_rOutMatrix.m[3][0] = 0.0f;
			_rOutMatrix.m[3][1] = 0.0f;
			_rOutMatrix.m[3][2] = 0.0f;
			_rOutMatrix.m[3][3] = 1.0f;
		}

		void Mtx_BuildRotationZ( Matrix44& _rOutMatrix, REAL _fAngleRadian )
		{
			REAL fSinAngle = Math::Sin( _fAngleRadian );
			REAL fCosAngle = Math::Cos( _fAngleRadian );

			_rOutMatrix.m[0][0] = fCosAngle;
			_rOutMatrix.m[0][1] = fSinAngle;
			_rOutMatrix.m[0][2] = 0.0f;
			_rOutMatrix.m[0][3] = 0.0f;

			_rOutMatrix.m[1][0] = -fSinAngle;
			_rOutMatrix.m[1][1] = fCosAngle;
			_rOutMatrix.m[1][2] = 0.0f;
			_rOutMatrix.m[1][3] = 0.0f;

			_rOutMatrix.m[2][0] = 0.0f;
			_rOutMatrix.m[2][1] = 0.0f;
			_rOutMatrix.m[2][2] = 1.0f;
			_rOutMatrix.m[2][3] = 0.0f;

			_rOutMatrix.m[3][0] = 0.0f;
			_rOutMatrix.m[3][1] = 0.0f;
			_rOutMatrix.m[3][2] = 0.0f;
			_rOutMatrix.m[3][3] = 1.0f;
		}
		void Mtx_BuildRotation( Matrix44& _rOutM, const Vec3& _rvAxis, REAL _fAngleRadian )
		{
			// XMVECTOR xmvAxis = XMVectorSet( vecAxis.GetX(), vecAxis.GetY(), vecAxis.GetZ(), 0.f );
			//
			// Matrix44 mtxRotation = XMMatrixRotationAxis( xmvAxis, Math::DegreeToRadian( fAngleDegree ) );
			//
			// matrix = XMMatrixMultiply( matrix, mtxRotation );

			float u = _rvAxis[0];
			float v = _rvAxis[1];
			float w = _rvAxis[2];

			float rcos = Math::Cos(_fAngleRadian);
			float rsin = Math::Sin(_fAngleRadian);

			//_rOutM.m[0][0] =      rcos + u*u*(1-rcos);
			//_rOutM.m[1][0] =  w * rsin + v*u*(1-rcos);
			//_rOutM.m[2][0] = -v * rsin + w*u*(1-rcos);

			//_rOutM.m[0][1] = -w * rsin + u*v*(1-rcos);
			//_rOutM.m[1][1] =      rcos + v*v*(1-rcos);
			//_rOutM.m[2][1] =  u * rsin + w*v*(1-rcos);

			//_rOutM.m[0][2] =  v * rsin + u*w*(1-rcos);
			//_rOutM.m[1][2] = -u * rsin + v*w*(1-rcos);
			//_rOutM.m[2][2] =      rcos + w*w*(1-rcos);

			_rOutM.m[0][0] =      rcos + u*u*(1-rcos);
			_rOutM.m[0][1] =  w * rsin + v*u*(1-rcos);
			_rOutM.m[0][2] = -v * rsin + w*u*(1-rcos);

			_rOutM.m[1][0] = -w * rsin + u*v*(1-rcos);
			_rOutM.m[1][1] =      rcos + v*v*(1-rcos);
			_rOutM.m[1][2] =  u * rsin + w*v*(1-rcos);

			_rOutM.m[2][0] =  v * rsin + u*w*(1-rcos);
			_rOutM.m[2][1] = -u * rsin + v*w*(1-rcos);
			_rOutM.m[2][2] =      rcos + w*w*(1-rcos);

			_rOutM.m[0][3] = 0.0f;
			_rOutM.m[1][3] = 0.0f;
			_rOutM.m[2][3] = 0.0f;

			_rOutM.m[3][0] = 0.0f;
			_rOutM.m[3][1] = 0.0f;
			_rOutM.m[3][2] = 0.0f;

			_rOutM.m[3][3] = 1.0f;

		}

				
		void Mtx_BuildFromSRT( Matrix44& _rMatrix, const Vec3 & _vTranslation, const Quaternion & _quatRotation, const Vec3 & _vScale )
		{
			// rotation part
			Matrix44 mtxRotation;
			_quatRotation.ToMatrix( mtxRotation );
			
			// scaling part
			Matrix44 mtxScale;
			Mtx_BuildScale( mtxScale, _vScale );
			Mtx_Concat( mtxScale, mtxRotation, _rMatrix );
			
			// simply set the translation coefficients
			_rMatrix._41 = _vTranslation[ 0 ];
			_rMatrix._42 = _vTranslation[ 1 ];
			_rMatrix._43 = _vTranslation[ 2 ];
		}
		
		void Mtx_BuildFromSRT( Matrix44& _rOutMatrix, const SRT& _rSRT )
		{
			Mtx_BuildFromSRT( _rOutMatrix, _rSRT.mvTranslate, _rSRT.mquatRotation, _rSRT.mvScale );
		}
		void Mtx_BuildPerspectiveRH( Matrix44& _rOutResult, REAL _fFovAngleY, REAL _fAspectRatio, REAL _fNearZ, REAL _fFarZ )
		{
			REAL    fSinFov;
			REAL	fCosFov;
			REAL	fHeight;
			REAL	fWidth; 

			//XMASSERT(!XMScalarNearEqual(_fFovAngleY, 0.0f, 0.00001f * 2.0f));
			//XMASSERT(!XMScalarNearEqual(_fAspectRatio, 0.0f, 0.00001f));
			//XMASSERT(!XMScalarNearEqual(FarZ, NearZ, 0.00001f));
 
			fSinFov = Math::Sin( 0.5f * _fFovAngleY );
			fCosFov = Math::Cos( 0.5f * _fFovAngleY );

			fHeight = fCosFov / fSinFov;
			fWidth = fHeight / _fAspectRatio;

			_rOutResult.m[0][0] = fWidth;
			_rOutResult.m[0][1] = 0.0f;
			_rOutResult.m[0][2] = 0.0f;
			_rOutResult.m[0][3] = 0.0f;

			_rOutResult.m[1][0] = 0.0f;
			_rOutResult.m[1][1] = fHeight ;
			_rOutResult.m[1][2] = 0.0f;
			_rOutResult.m[1][3] = 0.0f;

 			_rOutResult.m[2][0] = 0.0f;
			_rOutResult.m[2][1] = 0.0f;
			_rOutResult.m[2][2] = _fFarZ / (_fNearZ - _fFarZ);
			_rOutResult.m[2][3] = -1.0f;

 			_rOutResult.m[3][0] = 0.0f;
			_rOutResult.m[3][1] = 0.0f;
			_rOutResult.m[3][2] = _rOutResult.m[2][2] * _fNearZ ;
			_rOutResult.m[3][3] = 0.0f;
 		}

		void Mtx_BuildOrtographicRH( Matrix44& _rOutResult, REAL _fViewWidth, REAL _fViewHeight, REAL _fNearZ, REAL _fFarZ )
		{
			//MB_ASSERT( );

			_rOutResult.m[0][0] = 2.0f / _fViewWidth;
			_rOutResult.m[0][1] = 0.0f;
			_rOutResult.m[0][2] = 0.0f;
			_rOutResult.m[0][3] = 0.0f;

			_rOutResult.m[1][0] = 0.0f;
			_rOutResult.m[1][1] = 2.0f / _fViewHeight;
			_rOutResult.m[1][2] = 0.0f;
			_rOutResult.m[1][3] = 0.0f;

 			_rOutResult.m[2][0] = 0.0f;
			_rOutResult.m[2][1] = 0.0f;
			_rOutResult.m[2][2] = 1.0f / (_fNearZ - _fFarZ);
			_rOutResult.m[2][3] = 0.0f;

 			_rOutResult.m[3][0] = 0.0f;
			_rOutResult.m[3][1] = 0.0f;
			_rOutResult.m[3][2] = _rOutResult.m[2][2] * _fNearZ ;
			_rOutResult.m[3][3] = 1.0f;
		}


		

		void Mtx_BuildFromSRT( Matrix43& _rMatrix, const SRT& _rSRT )
		{
 		//	// rotation part
			//_quatRotation.ToMatrix( _rMatrix );
			//
			//// scaling part
			//Mtx_Scale( _rMatrix, _vScale );
			//
			//// simply set the translation coefficients
			//Mtx_Translation( _rMatrix, _vTranslation );
			///*m[12] = _translation.x;
			//m[13] = _translation.y;
			//m[14] = _translation.z;*/

			Matrix44 tmp;

			Mtx_BuildIdentity( tmp );

			tmp._11 = _rMatrix._11;
			tmp._12 = _rMatrix._12;
			tmp._13 = _rMatrix._13;

			tmp._21 = _rMatrix._21;
			tmp._22 = _rMatrix._22;
			tmp._23 = _rMatrix._23;
 
			tmp._31 = _rMatrix._31;
			tmp._32 = _rMatrix._32;
			tmp._33 = _rMatrix._33;
 
			tmp._41 = _rMatrix._41;
			tmp._42 = _rMatrix._42;
			tmp._43 = _rMatrix._43;

			Mtx_BuildFromSRT( tmp, _rSRT );

			_rMatrix._11 = tmp._11;
			_rMatrix._12 = tmp._12;
			_rMatrix._13 = tmp._13;

			_rMatrix._21 = tmp._21;
			_rMatrix._22 = tmp._22;
			_rMatrix._23 = tmp._23;
 
			_rMatrix._31 = tmp._31;
			_rMatrix._32 = tmp._32;
			_rMatrix._33 = tmp._33;
 
			_rMatrix._41 = tmp._41;
			_rMatrix._42 = tmp._42;
			_rMatrix._43 = tmp._43;
		}
		void Mtx_MultVector( const Matrix44& _rM, const Vec3& _rVec, Vec3& _rvOutResult )
		{
			//XMVECTOR xmvVec = XMVectorSet( vec.GetX(), vec.GetY(), vec.GetZ(), 0.f );
			//XMVECTOR xmvVecResult = XMVector3Transform( xmvVec, matrix );
			//
			//vecResult.SetX( XMVectorGetX( xmvVecResult ) );
			//vecResult.SetY( XMVectorGetY( xmvVecResult ) );
			//vecResult.SetZ( XMVectorGetZ( xmvVecResult ) );
			//

 			_rvOutResult[ 0 ] =	_rVec[ 0 ] * _rM.m[ 0 ][ 0 ] 
							+	_rVec[ 1 ] * _rM.m[ 1 ][ 0 ] 
							+	_rVec[ 2 ] * _rM.m[ 2 ][ 0 ] 
							+	/* 1* */	_rM.m[ 3 ][ 0 ];

			_rvOutResult[ 1 ] =	_rVec[ 0 ] * _rM.m[ 0 ][ 1 ] 
							+	_rVec[ 1 ] * _rM.m[ 1 ][ 1 ] 
							+	_rVec[ 2 ] * _rM.m[ 2 ][ 1 ] 
							+	/* 1* */	_rM.m[ 3 ][ 1 ];

			_rvOutResult[ 2 ] =	_rVec[ 0 ] * _rM.m[ 0 ][ 2 ] 
							+	_rVec[ 1 ] * _rM.m[ 1 ][ 2 ] 
							+	_rVec[ 2 ] * _rM.m[ 2 ][ 2 ] 
							+	/* 1* */	_rM.m[ 3 ][ 2 ];
		}

		void Mtx_MultVector( const Matrix33& _rMatrix, const Vec3& vec, Vec3& vecResult )
		{
 			vecResult[ 0 ] =	vec[ 0 ] * _rMatrix.m[ 0 ][ 0 ] 
							+	vec[ 1 ] * _rMatrix.m[ 1 ][ 0 ] 
							+	vec[ 2 ] * _rMatrix.m[ 2 ][ 0 ];

			vecResult[ 1 ] =	vec[ 0 ] * _rMatrix.m[ 0 ][ 1 ] 
							+	vec[ 1 ] * _rMatrix.m[ 1 ][ 1 ] 
							+	vec[ 2 ] * _rMatrix.m[ 2 ][ 1 ];

			vecResult[ 2 ] =	vec[ 0 ] * _rMatrix.m[ 0 ][ 2 ] 
							+	vec[ 1 ] * _rMatrix.m[ 1 ][ 2 ] 
							+	vec[ 2 ] * _rMatrix.m[ 2 ][ 2 ];
    	}
		
		void Mtx_Concat( const Matrix44& _rM1, const Matrix44& _rM2, Matrix44& _rOutMatrix )
		{
			MB_ASSERT( & _rOutMatrix != & _rM1, "Cannot multiply 2 matrices in place" );
			MB_ASSERT( & _rOutMatrix != & _rM2, "Cannot multiply 2 matrices in place" );

			
			// Cache the invariants in registers
			float x = _rM1.m[0][0];
			float y = _rM1.m[0][1];
			float z = _rM1.m[0][2];
			float w = _rM1.m[0][3];
			// Perform the operation on the first row
			_rOutMatrix.m[0][0] = (_rM2.m[0][0]*x)+(_rM2.m[1][0]*y)+(_rM2.m[2][0]*z)+(_rM2.m[3][0]*w);
			_rOutMatrix.m[0][1] = (_rM2.m[0][1]*x)+(_rM2.m[1][1]*y)+(_rM2.m[2][1]*z)+(_rM2.m[3][1]*w);
			_rOutMatrix.m[0][2] = (_rM2.m[0][2]*x)+(_rM2.m[1][2]*y)+(_rM2.m[2][2]*z)+(_rM2.m[3][2]*w);
			_rOutMatrix.m[0][3] = (_rM2.m[0][3]*x)+(_rM2.m[1][3]*y)+(_rM2.m[2][3]*z)+(_rM2.m[3][3]*w);
			// Repeat for all the other rows
			x = _rM1.m[1][0];
			y = _rM1.m[1][1];
			z = _rM1.m[1][2];
			w = _rM1.m[1][3];
			_rOutMatrix.m[1][0] = (_rM2.m[0][0]*x)+(_rM2.m[1][0]*y)+(_rM2.m[2][0]*z)+(_rM2.m[3][0]*w);
			_rOutMatrix.m[1][1] = (_rM2.m[0][1]*x)+(_rM2.m[1][1]*y)+(_rM2.m[2][1]*z)+(_rM2.m[3][1]*w);
			_rOutMatrix.m[1][2] = (_rM2.m[0][2]*x)+(_rM2.m[1][2]*y)+(_rM2.m[2][2]*z)+(_rM2.m[3][2]*w);
			_rOutMatrix.m[1][3] = (_rM2.m[0][3]*x)+(_rM2.m[1][3]*y)+(_rM2.m[2][3]*z)+(_rM2.m[3][3]*w);
			x = _rM1.m[2][0];
			y = _rM1.m[2][1];
			z = _rM1.m[2][2];
			w = _rM1.m[2][3];
			_rOutMatrix.m[2][0] = (_rM2.m[0][0]*x)+(_rM2.m[1][0]*y)+(_rM2.m[2][0]*z)+(_rM2.m[3][0]*w);
			_rOutMatrix.m[2][1] = (_rM2.m[0][1]*x)+(_rM2.m[1][1]*y)+(_rM2.m[2][1]*z)+(_rM2.m[3][1]*w);
			_rOutMatrix.m[2][2] = (_rM2.m[0][2]*x)+(_rM2.m[1][2]*y)+(_rM2.m[2][2]*z)+(_rM2.m[3][2]*w);
			_rOutMatrix.m[2][3] = (_rM2.m[0][3]*x)+(_rM2.m[1][3]*y)+(_rM2.m[2][3]*z)+(_rM2.m[3][3]*w);
			x = _rM1.m[3][0];
			y = _rM1.m[3][1];
			z = _rM1.m[3][2];
			w = _rM1.m[3][3];
			_rOutMatrix.m[3][0] = (_rM2.m[0][0]*x)+(_rM2.m[1][0]*y)+(_rM2.m[2][0]*z)+(_rM2.m[3][0]*w);
			_rOutMatrix.m[3][1] = (_rM2.m[0][1]*x)+(_rM2.m[1][1]*y)+(_rM2.m[2][1]*z)+(_rM2.m[3][1]*w);
			_rOutMatrix.m[3][2] = (_rM2.m[0][2]*x)+(_rM2.m[1][2]*y)+(_rM2.m[2][2]*z)+(_rM2.m[3][2]*w);
			_rOutMatrix.m[3][3] = (_rM2.m[0][3]*x)+(_rM2.m[1][3]*y)+(_rM2.m[2][3]*z)+(_rM2.m[3][3]*w);
			// NOTE: might only works if Matrix44 is XMMATRIX...
			//matrixResult = XMMatrixMultiply( matrixOne, matrixTwo ); 
		}
 
		void Mtx_Add( const Matrix44& matrixOne, const Matrix44& matrixTwo, Matrix44& matrixResult )
		{
			matrixResult.m[0][0] = matrixOne.m[0][0] + matrixTwo.m[0][0];
			matrixResult.m[0][1] = matrixOne.m[0][1] + matrixTwo.m[0][1];
			matrixResult.m[0][2] = matrixOne.m[0][2] + matrixTwo.m[0][2];
			matrixResult.m[0][3] = matrixOne.m[0][3] + matrixTwo.m[0][3];
 
			matrixResult.m[1][0] = matrixOne.m[1][0] + matrixTwo.m[1][0];
			matrixResult.m[1][1] = matrixOne.m[1][1] + matrixTwo.m[1][1];
			matrixResult.m[1][2] = matrixOne.m[1][2] + matrixTwo.m[1][2];
			matrixResult.m[1][3] = matrixOne.m[1][3] + matrixTwo.m[1][3];
 
			matrixResult.m[2][0] = matrixOne.m[2][0] + matrixTwo.m[2][0];
			matrixResult.m[2][1] = matrixOne.m[2][1] + matrixTwo.m[2][1];
			matrixResult.m[2][2] = matrixOne.m[2][2] + matrixTwo.m[2][2];
			matrixResult.m[2][3] = matrixOne.m[2][3] + matrixTwo.m[2][3];
 
			matrixResult.m[3][0] = matrixOne.m[3][0] + matrixTwo.m[3][0];
			matrixResult.m[3][1] = matrixOne.m[3][1] + matrixTwo.m[3][1];
			matrixResult.m[3][2] = matrixOne.m[3][2] + matrixTwo.m[3][2];
  			matrixResult.m[3][3] = matrixOne.m[3][3] + matrixTwo.m[3][3];
  		}


		void Mtx_BuildNull( Matrix43& _rOutMatrix )
		{
 			_rOutMatrix.m[0][0] = 0.0f;
			_rOutMatrix.m[0][1] = 0.0f;
			_rOutMatrix.m[0][2] = 0.0f;
 
			_rOutMatrix.m[1][0] = 0.0f;
			_rOutMatrix.m[1][1] = 0.0f;
			_rOutMatrix.m[1][2] = 0.0f;
 
			_rOutMatrix.m[2][0] = 0.0f;
			_rOutMatrix.m[2][1] = 0.0f;
			_rOutMatrix.m[2][2] = 0.0f;
 
			_rOutMatrix.m[3][0] = 0.0f;
			_rOutMatrix.m[3][1] = 0.0f;
			_rOutMatrix.m[3][2] = 0.0f;
 		}

		
		void Mtx_Transpose( const Matrix44& _rMatrix, Matrix44& _rOutResult )
		{			
			MB_ASSERT( & _rMatrix != & _rOutResult, "Cannot transpose in place" );

			//_rOutResult = XMMatrixTranspose( _rMatrix );

			// diagonal stays the same
			_rOutResult._11 = _rMatrix._11;	
			_rOutResult._12 = _rMatrix._21;
			_rOutResult._13 = _rMatrix._31;
			_rOutResult._14 = _rMatrix._41;
			
			_rOutResult._21 = _rMatrix._12;
			_rOutResult._22 = _rMatrix._22;
			_rOutResult._23 = _rMatrix._32;
			_rOutResult._24 = _rMatrix._42;

			_rOutResult._31 = _rMatrix._13;
			_rOutResult._32 = _rMatrix._23;
			_rOutResult._33 = _rMatrix._33;
			_rOutResult._34 = _rMatrix._43;

			_rOutResult._41 = _rMatrix._14;
			_rOutResult._42 = _rMatrix._24;
			_rOutResult._43 = _rMatrix._34;
			_rOutResult._44 = _rMatrix._44;
		}

		void Mtx_Inverse( const Matrix44& _rMatrix, Matrix44& _rOutResult , Vec3& _rvOutDeterminant )
		{
			//XMVECTOR xmVec;
			//_rOutResult = XMMatrixInverse( &xmVec, _rMatrix );			
			//_rvOutDeterminant.SetX( XMVectorGetX( xmVec ) );
			//_rvOutDeterminant.SetY( XMVectorGetY( xmVec ) );
			//_rvOutDeterminant.SetZ( XMVectorGetZ( xmVec ) );

 			float inv[16], det;

			inv[0] = _rMatrix.ma[5]  * _rMatrix.ma[10] * _rMatrix.ma[15] - 
					 _rMatrix.ma[5]  * _rMatrix.ma[11] * _rMatrix.ma[14] - 
					 _rMatrix.ma[9]  * _rMatrix.ma[6]  * _rMatrix.ma[15] + 
					 _rMatrix.ma[9]  * _rMatrix.ma[7]  * _rMatrix.ma[14] +
					 _rMatrix.ma[13] * _rMatrix.ma[6]  * _rMatrix.ma[11] - 
					 _rMatrix.ma[13] * _rMatrix.ma[7]  * _rMatrix.ma[10];

			inv[4] = -_rMatrix.ma[4]  * _rMatrix.ma[10] * _rMatrix.ma[15] + 
					  _rMatrix.ma[4]  * _rMatrix.ma[11] * _rMatrix.ma[14] + 
					  _rMatrix.ma[8]  * _rMatrix.ma[6]  * _rMatrix.ma[15] - 
					  _rMatrix.ma[8]  * _rMatrix.ma[7]  * _rMatrix.ma[14] - 
					  _rMatrix.ma[12] * _rMatrix.ma[6]  * _rMatrix.ma[11] + 
					  _rMatrix.ma[12] * _rMatrix.ma[7]  * _rMatrix.ma[10];

			inv[8] = _rMatrix.ma[4]  * _rMatrix.ma[9] * _rMatrix.ma[15] - 
					 _rMatrix.ma[4]  * _rMatrix.ma[11] * _rMatrix.ma[13] - 
					 _rMatrix.ma[8]  * _rMatrix.ma[5] * _rMatrix.ma[15] + 
					 _rMatrix.ma[8]  * _rMatrix.ma[7] * _rMatrix.ma[13] + 
					 _rMatrix.ma[12] * _rMatrix.ma[5] * _rMatrix.ma[11] - 
					 _rMatrix.ma[12] * _rMatrix.ma[7] * _rMatrix.ma[9];

			inv[12] = -_rMatrix.ma[4]  * _rMatrix.ma[9] * _rMatrix.ma[14] + 
					   _rMatrix.ma[4]  * _rMatrix.ma[10] * _rMatrix.ma[13] +
					   _rMatrix.ma[8]  * _rMatrix.ma[5] * _rMatrix.ma[14] - 
					   _rMatrix.ma[8]  * _rMatrix.ma[6] * _rMatrix.ma[13] - 
					   _rMatrix.ma[12] * _rMatrix.ma[5] * _rMatrix.ma[10] + 
					   _rMatrix.ma[12] * _rMatrix.ma[6] * _rMatrix.ma[9];

			inv[1] = -_rMatrix.ma[1]  * _rMatrix.ma[10] * _rMatrix.ma[15] + 
					  _rMatrix.ma[1]  * _rMatrix.ma[11] * _rMatrix.ma[14] + 
					  _rMatrix.ma[9]  * _rMatrix.ma[2] * _rMatrix.ma[15] - 
					  _rMatrix.ma[9]  * _rMatrix.ma[3] * _rMatrix.ma[14] - 
					  _rMatrix.ma[13] * _rMatrix.ma[2] * _rMatrix.ma[11] + 
					  _rMatrix.ma[13] * _rMatrix.ma[3] * _rMatrix.ma[10];

			inv[5] = _rMatrix.ma[0]  * _rMatrix.ma[10] * _rMatrix.ma[15] - 
					 _rMatrix.ma[0]  * _rMatrix.ma[11] * _rMatrix.ma[14] - 
					 _rMatrix.ma[8]  * _rMatrix.ma[2] * _rMatrix.ma[15] + 
					 _rMatrix.ma[8]  * _rMatrix.ma[3] * _rMatrix.ma[14] + 
					 _rMatrix.ma[12] * _rMatrix.ma[2] * _rMatrix.ma[11] - 
					 _rMatrix.ma[12] * _rMatrix.ma[3] * _rMatrix.ma[10];

			inv[9] = -_rMatrix.ma[0]  * _rMatrix.ma[9] * _rMatrix.ma[15] + 
					  _rMatrix.ma[0]  * _rMatrix.ma[11] * _rMatrix.ma[13] + 
					  _rMatrix.ma[8]  * _rMatrix.ma[1] * _rMatrix.ma[15] - 
					  _rMatrix.ma[8]  * _rMatrix.ma[3] * _rMatrix.ma[13] - 
					  _rMatrix.ma[12] * _rMatrix.ma[1] * _rMatrix.ma[11] + 
					  _rMatrix.ma[12] * _rMatrix.ma[3] * _rMatrix.ma[9];

			inv[13] = _rMatrix.ma[0]  * _rMatrix.ma[9] * _rMatrix.ma[14] - 
					  _rMatrix.ma[0]  * _rMatrix.ma[10] * _rMatrix.ma[13] - 
					  _rMatrix.ma[8]  * _rMatrix.ma[1] * _rMatrix.ma[14] + 
					  _rMatrix.ma[8]  * _rMatrix.ma[2] * _rMatrix.ma[13] + 
					  _rMatrix.ma[12] * _rMatrix.ma[1] * _rMatrix.ma[10] - 
					  _rMatrix.ma[12] * _rMatrix.ma[2] * _rMatrix.ma[9];

			inv[2] = _rMatrix.ma[1]  * _rMatrix.ma[6] * _rMatrix.ma[15] - 
					 _rMatrix.ma[1]  * _rMatrix.ma[7] * _rMatrix.ma[14] - 
					 _rMatrix.ma[5]  * _rMatrix.ma[2] * _rMatrix.ma[15] + 
					 _rMatrix.ma[5]  * _rMatrix.ma[3] * _rMatrix.ma[14] + 
					 _rMatrix.ma[13] * _rMatrix.ma[2] * _rMatrix.ma[7] - 
					 _rMatrix.ma[13] * _rMatrix.ma[3] * _rMatrix.ma[6];

			inv[6] = -_rMatrix.ma[0]  * _rMatrix.ma[6] * _rMatrix.ma[15] + 
					  _rMatrix.ma[0]  * _rMatrix.ma[7] * _rMatrix.ma[14] + 
					  _rMatrix.ma[4]  * _rMatrix.ma[2] * _rMatrix.ma[15] - 
					  _rMatrix.ma[4]  * _rMatrix.ma[3] * _rMatrix.ma[14] - 
					  _rMatrix.ma[12] * _rMatrix.ma[2] * _rMatrix.ma[7] + 
					  _rMatrix.ma[12] * _rMatrix.ma[3] * _rMatrix.ma[6];

			inv[10] = _rMatrix.ma[0]  * _rMatrix.ma[5] * _rMatrix.ma[15] - 
					  _rMatrix.ma[0]  * _rMatrix.ma[7] * _rMatrix.ma[13] - 
					  _rMatrix.ma[4]  * _rMatrix.ma[1] * _rMatrix.ma[15] + 
					  _rMatrix.ma[4]  * _rMatrix.ma[3] * _rMatrix.ma[13] + 
					  _rMatrix.ma[12] * _rMatrix.ma[1] * _rMatrix.ma[7] - 
					  _rMatrix.ma[12] * _rMatrix.ma[3] * _rMatrix.ma[5];

			inv[14] = -_rMatrix.ma[0]  * _rMatrix.ma[5] * _rMatrix.ma[14] + 
					   _rMatrix.ma[0]  * _rMatrix.ma[6] * _rMatrix.ma[13] + 
					   _rMatrix.ma[4]  * _rMatrix.ma[1] * _rMatrix.ma[14] - 
					   _rMatrix.ma[4]  * _rMatrix.ma[2] * _rMatrix.ma[13] - 
					   _rMatrix.ma[12] * _rMatrix.ma[1] * _rMatrix.ma[6] + 
					   _rMatrix.ma[12] * _rMatrix.ma[2] * _rMatrix.ma[5];

			inv[3] = -_rMatrix.ma[1] * _rMatrix.ma[6] * _rMatrix.ma[11] + 
					  _rMatrix.ma[1] * _rMatrix.ma[7] * _rMatrix.ma[10] + 
					  _rMatrix.ma[5] * _rMatrix.ma[2] * _rMatrix.ma[11] - 
					  _rMatrix.ma[5] * _rMatrix.ma[3] * _rMatrix.ma[10] - 
					  _rMatrix.ma[9] * _rMatrix.ma[2] * _rMatrix.ma[7] + 
					  _rMatrix.ma[9] * _rMatrix.ma[3] * _rMatrix.ma[6];

			inv[7] = _rMatrix.ma[0] * _rMatrix.ma[6] * _rMatrix.ma[11] - 
					 _rMatrix.ma[0] * _rMatrix.ma[7] * _rMatrix.ma[10] - 
					 _rMatrix.ma[4] * _rMatrix.ma[2] * _rMatrix.ma[11] + 
					 _rMatrix.ma[4] * _rMatrix.ma[3] * _rMatrix.ma[10] + 
					 _rMatrix.ma[8] * _rMatrix.ma[2] * _rMatrix.ma[7] - 
					 _rMatrix.ma[8] * _rMatrix.ma[3] * _rMatrix.ma[6];

			inv[11] = -_rMatrix.ma[0] * _rMatrix.ma[5] * _rMatrix.ma[11] + 
					   _rMatrix.ma[0] * _rMatrix.ma[7] * _rMatrix.ma[9] + 
					   _rMatrix.ma[4] * _rMatrix.ma[1] * _rMatrix.ma[11] - 
					   _rMatrix.ma[4] * _rMatrix.ma[3] * _rMatrix.ma[9] - 
					   _rMatrix.ma[8] * _rMatrix.ma[1] * _rMatrix.ma[7] + 
					   _rMatrix.ma[8] * _rMatrix.ma[3] * _rMatrix.ma[5];

			inv[15] = _rMatrix.ma[0] * _rMatrix.ma[5] * _rMatrix.ma[10] - 
					  _rMatrix.ma[0] * _rMatrix.ma[6] * _rMatrix.ma[9] - 
					  _rMatrix.ma[4] * _rMatrix.ma[1] * _rMatrix.ma[10] + 
					  _rMatrix.ma[4] * _rMatrix.ma[2] * _rMatrix.ma[9] + 
					  _rMatrix.ma[8] * _rMatrix.ma[1] * _rMatrix.ma[6] - 
					  _rMatrix.ma[8] * _rMatrix.ma[2] * _rMatrix.ma[5];

			det = _rMatrix.ma[0] * inv[0] + _rMatrix.ma[1] * inv[4] + _rMatrix.ma[2] * inv[8] + _rMatrix.ma[3] * inv[12];

			_rvOutDeterminant[ 0 ] = det;
			_rvOutDeterminant[ 1 ] = det;
			_rvOutDeterminant[ 2 ] = det;

			if( det != 0 )
			{
				det = 1.0f / det;

				for( int i = 0; i < 16; i++)
					_rOutResult.ma[i] = inv[i] * det;

			}
		}
		
	}
}


//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//	Implementations :	QUATERNION 
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
namespace Melting
{
	namespace Math
	{
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion::Quaternion( void )
		:	mfRealPart		()
		,	mvImaginaryPart	()
		{
			// !!! KEEP IT UNINITIALIZED FOR PERFORMANCE REASON !!!!
		}
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion::Quaternion( REAL _fRealPart, REAL _fImaginaryPartX, REAL _fImaginaryPartY, REAL _fImaginaryPartZ )
		:	mfRealPart		( _fRealPart )
		,	mvImaginaryPart	( _fImaginaryPartX, _fImaginaryPartY, _fImaginaryPartZ )
		{
		}
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion::Quaternion( const Quaternion & _rQuat )
		:	mfRealPart		( _rQuat.mfRealPart )
		,	mvImaginaryPart	( _rQuat.mvImaginaryPart )
		{ 
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion::Quaternion( Radian rAngle, const Vec3 & _rvAxis )
		:	mfRealPart		()
		,	mvImaginaryPart	()
		{
			REAL fSin( Math::Sin( rAngle * .5f ) );

			mvImaginaryPart.SetX( _rvAxis.GetX() * fSin );
			mvImaginaryPart.SetY( _rvAxis.GetY() * fSin );
			mvImaginaryPart.SetZ( _rvAxis.GetZ() * fSin );
			mfRealPart = Math::Cos( rAngle * .5f );

			if( GetLength() > 0.f )
			{
				Normalize();
			}
		}
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 		
		Quaternion::Quaternion( REAL _angleX, REAL _angleY, REAL _angleZ)
		{
			FromAngles( _angleX, _angleY, _angleZ );
		}
		
		//---------------------------------------------------------------------------------------------- 
		// Construct a Quaternion from the given Euler Angles
		//---------------------------------------------------------------------------------------------- 
		void Quaternion::FromAngles	( const Vec3 & _rvAngles )
		{
			FromAngles( _rvAngles[ 0 ], _rvAngles[ 1 ], _rvAngles[ 2 ] );
		}

		//---------------------------------------------------------------------------------------------- 
		// Construct a Quaternion from the given Euler Angles
		//---------------------------------------------------------------------------------------------- 
		void Quaternion::FromAngles	( REAL _angleX, REAL _angleY, REAL _angleZ )
		{
			REAL angle;
			REAL sr, sp, sy, cr, cp, cy;
			
			angle 	= _angleZ * REAL( 0.5f );
			sy 		= Math::Sin(angle);
			cy 		= Math::Cos(angle);
		
			angle 	= _angleY * REAL( 0.5f );
			sp 		= Math::Sin(angle);
			cp 		= Math::Cos(angle);
		
			angle 	= _angleX * REAL( 0.5f );
			sr 		= Math::Sin(angle);
			cr 		= Math::Cos(angle);
		
			REAL crcp = cr*cp;
			REAL srsp = sr*sp;
		
			mvImaginaryPart[0] = ( sr*cp*cy	- cr*sp*sy );
			mvImaginaryPart[1] = ( cr*sp*cy	+ sr*cp*sy );
			mvImaginaryPart[2] = ( crcp*sy	- srsp*cy );
			mfRealPart = ( crcp*cy	+ srsp*sy ); 
		} 

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		void Quaternion::FromDirections	( const Vec3 & _rvOldDirection, const Vec3 & _rvNewDirection )
		{
			// If V1 and V2 are not parallel, the axis of rotation is the unit-length
			// vector U = Cross(V1,V2)/Length(Cross(V1,V2)).  The angle of rotation,
			// A, is the angle between V1 and V2.  The quaternion for the rotation is
			// q = cos(A/2) + sin(A/2)*(ux*i+uy*j+uz*k) where U = (ux,uy,uz).
			//
			// (1) Rather than extract A = acos(Dot(V1,V2)), multiply by 1/2, then
			//     compute sin(A/2) and cos(A/2), we reduce the computational costs by
			//     computing the bisector B = (V1+V2)/Length(V1+V2), so cos(A/2) =
			//     Dot(V1,B).
			//
			// (2) The rotation axis is U = Cross(V1,B)/Length(Cross(V1,B)), but
			//     Length(Cross(V1,B)) = Length(V1)*Length(B)*sin(A/2) = sin(A/2), in
			//     which case sin(A/2)*(ux*i+uy*j+uz*k) = (cx*i+cy*j+cz*k) where
			//     C = Cross(V1,B).
			//
			// If V1 = V2, then B = V1, cos(A/2) = 1, and U = (0,0,0).  If V1 = -V2,
			// then B = 0.  This can happen even if V1 is approximately -V2 using
			// floating point arithmetic, since Vector3::Normalize checks for
			// closeness to zero and returns the zero vector accordingly.  The test
			// for exactly zero is usually not recommend for floating point
			// arithmetic, but the implementation of Vector3::Normalize guarantees
			// the comparison is robust.  In this case, the A = pi and any axis
			// perpendicular to V1 may be used as the rotation axis.

			Vec3 vBisector = _rvOldDirection + _rvNewDirection;
			vBisector.Normalize();

			REAL fCosHalfAngle = _rvNewDirection.DotP( vBisector );
			Vec3 vCross;

			mfRealPart = fCosHalfAngle;

			if ( fCosHalfAngle != REAL( 0.0f ) )
			{
				mvImaginaryPart = _rvNewDirection.CrossP( vBisector );
			}
			else
			{
				REAL fInvLength;

				if ( Math::Abs(_rvNewDirection[0]) >= Math::Abs(_rvNewDirection[1]))
				{
					// V1.x or V1.z is the largest magnitude component
					fInvLength = Math::InverseSquareRoot	( _rvNewDirection[0] * _rvNewDirection[0] 
																+ _rvNewDirection[2] * _rvNewDirection[2]);

					mvImaginaryPart[0] = - _rvNewDirection[2] * fInvLength;
					mvImaginaryPart[1] = REAL( 0.0f );
					mvImaginaryPart[2] = + _rvNewDirection[0] * fInvLength;
				}
				else
				{
					// V1.y or V1.z is the largest magnitude component
					fInvLength = Math::InverseSquareRoot	( _rvNewDirection[1] * _rvNewDirection[1] 
																+ _rvNewDirection[2] * _rvNewDirection[2]);

					mvImaginaryPart[0] = REAL( 0.0f );
					mvImaginaryPart[1] = + _rvNewDirection[2] * fInvLength;
					mvImaginaryPart[2] = - _rvNewDirection[1] * fInvLength;
				}
			} 
		} 
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		void Quaternion::FromMatrix( const Matrix33& _rMatrix )
		{			
			//Calculate the trace of the matrix T from the equation:
			//           2     2     2
			//T = 4 - 4x  - 4y  - 4z
			//            2    2    2
			//   = 4( 1 -x  - y  - z )
			//   = 1 + mat[0] + mat[5] + mat[10]
			float T = _rMatrix.ma[0] + _rMatrix.ma[4] + _rMatrix.ma[8];

			//If the trace of the matrix is greater than zero, then
			//perform an "instant" calculation.
			//Important note wrt. rouning errors:
			// Test if ( T > 0.00000001 ) to avoid large distortions!

			float S, X, Y, Z, W;

			if( T > 0.0f ) 
			{
				S = sqrt( 1.0f + T ) * 2.0f;
				X = ( _rMatrix.ma[7] - _rMatrix.ma[5] ) / S;
				Y = ( _rMatrix.ma[2] - _rMatrix.ma[6] ) / S;
				Z = ( _rMatrix.ma[3] - _rMatrix.ma[1] ) / S;
				W = 0.25f * S;
			}
			else
			//If the trace of the matrix is equal to zero then identify
			//which major diagonal element has the greatest value.
			//Depending on this, calculate the following:
			if ( _rMatrix.ma[0] > _rMatrix.ma[4] && _rMatrix.ma[0] > _rMatrix.ma[8] )  // Column 0: 
			{	
				S  = sqrt( 1.0f + _rMatrix.ma[0] - _rMatrix.ma[4] - _rMatrix.ma[8] ) * 2.0f;
				X = 0.25f * S;
				Y = (_rMatrix.ma[3] + _rMatrix.ma[1] ) / S;
				Z = (_rMatrix.ma[2] + _rMatrix.ma[6] ) / S;
				W = (_rMatrix.ma[7] - _rMatrix.ma[5] ) / S;
			} 
			else 
			if ( _rMatrix.ma[4] > _rMatrix.ma[8] ) // Column 1: 
			{			
				S  = sqrt( 1.0f + _rMatrix.ma[4] - _rMatrix.ma[0] - _rMatrix.ma[8] ) * 2.0f;
				X = (_rMatrix.ma[3] + _rMatrix.ma[1] ) / S;
				Y = 0.25f * S;
				Z = (_rMatrix.ma[7] + _rMatrix.ma[5] ) / S;
				W = (_rMatrix.ma[2] - _rMatrix.ma[6] ) / S;
			} 
			else // Column 2:	
			{						
				S  = sqrt( 1.0f + _rMatrix.ma[8] - _rMatrix.ma[0] - _rMatrix.ma[4] ) * 2.0f;
				X = (_rMatrix.ma[2] + _rMatrix.ma[6] ) / S;
				Y = (_rMatrix.ma[7] + _rMatrix.ma[5] ) / S;
				Z = 0.25f * S;
				W = (_rMatrix.ma[3] - _rMatrix.ma[1] ) / S;
			}

			 //The quaternion is then defined as:
			 //  Q = | X Y Z W |
			this->mvImaginaryPart[0]	= X;
			this->mvImaginaryPart[1]	= Y;
			this->mvImaginaryPart[2]	= Z;
			this->mfRealPart			= W;
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Vec3 Quaternion::Rotate	( const Vec3 & _rvOldDirection ) const
		{
//			this->Normalize();
			Quaternion qPosInitial;
			qPosInitial.FromAngles( _rvOldDirection );
			Quaternion qPosFinal = (*this)* qPosInitial * this->GetInverse();

			return qPosFinal.GetImaginaryPart();
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 		
	    Vec3 Quaternion::operator * ( const Vec3& _rvVec ) const
		{ 
			// from nVidia SDK implementation
			Vec3 uv, uuv;
			Vec3 qvec( mvImaginaryPart );
			uv 	= qvec.CrossP(_rvVec);
			uuv = qvec.CrossP(uv);
			uv 	*= REAL( 2.f ) * mfRealPart;
			uuv *= REAL( 2.f );

			return _rvVec + uv + uuv;
	    }  
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		void Quaternion::ToMatrix( Matrix44 & _rOutMatrix ) const
		{ 
			// each tmp var will be used twice !
			REAL fXX = GetX() * GetX();
			REAL fXY = GetX() * GetY();
			REAL fXZ = GetX() * GetZ();
			REAL fXW = GetX() * GetW();

			REAL fYY = GetY() * GetY();
			REAL fYZ = GetY() * GetZ();
			REAL fYW = GetY() * GetW();

			REAL fZZ = GetZ() * GetZ();
			REAL fZW = GetZ() * GetW();

			_rOutMatrix.m[ 0 ][ 0 ] = REAL ( 1.f ) - REAL ( 2.f ) * ( fYY + fZZ );
			_rOutMatrix.m[ 0 ][ 1 ] = REAL ( 2.f ) * ( fXY + fZW );
			_rOutMatrix.m[ 0 ][ 2 ] = REAL ( 2.f ) * ( fXZ - fYW );
			_rOutMatrix.m[ 0 ][ 3 ] = REAL ( 0.f ) ;

			_rOutMatrix.m[ 1 ][ 0 ] = REAL ( 2.f ) * ( fXY - fZW );
			_rOutMatrix.m[ 1 ][ 1 ] = REAL ( 1.f ) - REAL ( 2.f ) * ( fXX + fZZ );
			_rOutMatrix.m[ 1 ][ 2 ] = REAL ( 2.f ) * ( fYZ + fXW );
			_rOutMatrix.m[ 1 ][ 3 ] = REAL ( 0.f ) ;

			_rOutMatrix.m[ 2 ][ 0 ] = REAL ( 2.f ) * ( fXZ + fYW );
			_rOutMatrix.m[ 2 ][ 1 ] = REAL ( 2.f ) * ( fYZ - fXW );
			_rOutMatrix.m[ 2 ][ 2 ] = REAL ( 1.f ) - REAL ( 2.f ) * ( fXX + fYY );
			_rOutMatrix.m[ 2 ][ 3 ] = REAL ( 0.f ) ;

			_rOutMatrix.m[ 3 ][ 0 ] = REAL ( 0.f ) ;
			_rOutMatrix.m[ 3 ][ 1 ] = REAL ( 0.f ) ;
			_rOutMatrix.m[ 3 ][ 2 ] = REAL ( 0.f ) ;
			_rOutMatrix.m[ 3 ][ 3 ] = REAL ( 1.f ) ;  

			//// D.R.Y. !
			//ToMatrix( (Matrix43&) _rOutMatrix );

			//_rOutMatrix.m[ 0 ][ 3 ] = REAL ( 0.f ) ;
			//_rOutMatrix.m[ 1 ][ 3 ] = REAL ( 0.f ) ;
			//_rOutMatrix.m[ 2 ][ 3 ] = REAL ( 0.f ) ;
			//_rOutMatrix.m[ 3 ][ 3 ] = REAL ( 1.f ) ;  

		}
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		void Quaternion::ToMatrix( Matrix43 & _rOutMatrix ) const
		{ 
			// each tmp var will be used twice !
			REAL fXX = GetX() * GetX();
			REAL fXY = GetX() * GetY();
			REAL fXZ = GetX() * GetZ();
			REAL fXW = GetX() * GetW();

			REAL fYY = GetY() * GetY();
			REAL fYZ = GetY() * GetZ();
			REAL fYW = GetY() * GetW();

			REAL fZZ = GetZ() * GetZ();
			REAL fZW = GetZ() * GetW();

			_rOutMatrix.m[ 0 ][ 0 ] = REAL ( 1.f ) - REAL ( 2.f ) * ( fYY + fZZ );
			_rOutMatrix.m[ 0 ][ 1 ] = REAL ( 2.f ) * ( fXY + fZW );
			_rOutMatrix.m[ 0 ][ 2 ] = REAL ( 2.f ) * ( fXZ - fYW );

			_rOutMatrix.m[ 1 ][ 0 ] = REAL ( 2.f ) * ( fXY - fZW );
			_rOutMatrix.m[ 1 ][ 1 ] = REAL ( 1.f ) - REAL ( 2.f ) * ( fXX + fZZ );
			_rOutMatrix.m[ 1 ][ 2 ] = REAL ( 2.f ) * ( fYZ + fXW );

			_rOutMatrix.m[ 2 ][ 0 ] = REAL ( 2.f ) * ( fXZ + fYW );
			_rOutMatrix.m[ 2 ][ 1 ] = REAL ( 2.f ) * ( fYZ - fXW );
			_rOutMatrix.m[ 2 ][ 2 ] = REAL ( 1.f ) - REAL ( 2.f ) * ( fXX + fYY );

			_rOutMatrix.m[ 3 ][ 0 ] = REAL ( 0.f ) ;
			_rOutMatrix.m[ 3 ][ 1 ] = REAL ( 0.f ) ;
			_rOutMatrix.m[ 3 ][ 2 ] = REAL ( 0.f ) ;

		}
		//----------------------------------------------------------------------------------------------  
		//---------------------------------------------------------------------------------------------- 
		REAL Quaternion::GetRealPart	( void ) const
		{
			return mfRealPart;
		}

		//----------------------------------------------------------------------------------------------  
		//---------------------------------------------------------------------------------------------- 
		Vec3 Quaternion::GetImaginaryPart( void ) const
		{
			return mvImaginaryPart;
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		REAL Quaternion::GetX( void ) const	
		{ 
			return mvImaginaryPart.GetX(); 
		}
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		REAL Quaternion::GetY( void ) const	
		{ 
			return mvImaginaryPart.GetY(); 
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		REAL Quaternion::GetZ( void ) const	
		{ 
			return mvImaginaryPart.GetZ(); 
		}
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		REAL Quaternion::GetW( void ) const	
		{ 
			return mfRealPart; 
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion Quaternion::GetInverse( void ) const
		{
			return Quaternion(mfRealPart, - mvImaginaryPart[0], - mvImaginaryPart[1], - mvImaginaryPart[2] );
		} 
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		void Quaternion::Inverse( const Quaternion& _rBase, Quaternion& _rInverse )
		{			
			_rInverse.mfRealPart			= _rBase.mfRealPart;
			_rInverse.mvImaginaryPart[0]	= - _rBase.mvImaginaryPart[0];
			_rInverse.mvImaginaryPart[1]	= - _rBase.mvImaginaryPart[1];
			_rInverse.mvImaginaryPart[2]	= - _rBase.mvImaginaryPart[2];

			
			//REAL fLen = _rBase.GetLength();
			//
			//_rInverse.mfRealPart			= _rBase.mfRealPart / fLen;
			//_rInverse.mvImaginaryPart[0]	= - _rBase.mvImaginaryPart[0] / fLen;
			//_rInverse.mvImaginaryPart[1]	= - _rBase.mvImaginaryPart[1] / fLen;
			//_rInverse.mvImaginaryPart[2]	= - _rBase.mvImaginaryPart[2] / fLen;
			
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		bool Quaternion::operator==( const Quaternion & _rOther )
		{
//			return	mfRealPart == _rOther.mfRealPart
			return	Math::Abs( mfRealPart - _rOther.mfRealPart) < 0.01
				&&	mvImaginaryPart == _rOther.mvImaginaryPart;
		}
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		bool Quaternion::operator!=( const Quaternion & _rOther )
		{
			return !( (*this) == _rOther );
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion Quaternion::operator+(  const Quaternion & _rQuat )	const
		{
		   return Quaternion( mfRealPart + _rQuat.mfRealPart
							, mvImaginaryPart[0] + _rQuat.mvImaginaryPart[0]
							, mvImaginaryPart[1] + _rQuat.mvImaginaryPart[1]
							, mvImaginaryPart[2] + _rQuat.mvImaginaryPart[2]);
		}
 
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion	Quaternion::operator-(  void )	const
		{  
		   return Quaternion( - mfRealPart, - mvImaginaryPart[0], - mvImaginaryPart[1], - mvImaginaryPart[2] );
		}
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion	Quaternion::operator-(  const Quaternion & _rQuat )	const
		{  
		   return Quaternion( mfRealPart - _rQuat.mfRealPart
							, mvImaginaryPart[0] - _rQuat.mvImaginaryPart[0]
							, mvImaginaryPart[1] - _rQuat.mvImaginaryPart[1]
							, mvImaginaryPart[2] - _rQuat.mvImaginaryPart[2]);
		}

		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion 	Quaternion::operator/( REAL _fScalar ) const
		{
			return Quaternion	( mfRealPart/_fScalar
								, mvImaginaryPart[0]/_fScalar
								, mvImaginaryPart[1]/_fScalar
								, mvImaginaryPart[2]/_fScalar);
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion 	Quaternion::operator*( REAL _fScalar )	const
		{
			return Quaternion	( mfRealPart * _fScalar
								, mvImaginaryPart[0] * _fScalar
								, mvImaginaryPart[1] * _fScalar
								, mvImaginaryPart[2] * _fScalar);
		}	
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 		
		Quaternion & Quaternion::operator+=( const Quaternion& _rQuat )
		{
			mvImaginaryPart[0]	+= _rQuat.mvImaginaryPart[0];
			mvImaginaryPart[1]	+= _rQuat.mvImaginaryPart[1];
			mvImaginaryPart[2]	+= _rQuat.mvImaginaryPart[2];
			mfRealPart			+= _rQuat.mfRealPart;
		
			return *this;
		}
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 		
		Quaternion & Quaternion::operator-=( const Quaternion& _rQuat )
		{
			mvImaginaryPart[0]	-= _rQuat.mvImaginaryPart[0];
			mvImaginaryPart[1]	-= _rQuat.mvImaginaryPart[1];
			mvImaginaryPart[2]	-= _rQuat.mvImaginaryPart[2];
			mfRealPart			-= _rQuat.mfRealPart;
			
			return *this;
		}  

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 		
		Quaternion Quaternion::operator*( const Quaternion &_rQuat ) const
		{
			   Quaternion ret;
		
			   ret.mvImaginaryPart[0]	= mfRealPart * _rQuat.mvImaginaryPart[0] + mvImaginaryPart[0] * _rQuat.mfRealPart + mvImaginaryPart[1] * _rQuat.mvImaginaryPart[2] - mvImaginaryPart[2] * _rQuat.mvImaginaryPart[1];
			   ret.mvImaginaryPart[1]	= mfRealPart * _rQuat.mvImaginaryPart[1] + mvImaginaryPart[1] * _rQuat.mfRealPart + mvImaginaryPart[2] * _rQuat.mvImaginaryPart[0] - mvImaginaryPart[0] * _rQuat.mvImaginaryPart[2];
			   ret.mvImaginaryPart[2]	= mfRealPart * _rQuat.mvImaginaryPart[2] + mvImaginaryPart[2] * _rQuat.mfRealPart + mvImaginaryPart[0] * _rQuat.mvImaginaryPart[1] - mvImaginaryPart[1] * _rQuat.mvImaginaryPart[0];
			   ret.mfRealPart			= mfRealPart * _rQuat.mfRealPart - mvImaginaryPart[0] * _rQuat.mvImaginaryPart[0] - mvImaginaryPart[1] * _rQuat.mvImaginaryPart[1] - mvImaginaryPart[2] * _rQuat.mvImaginaryPart[2];

			   return ret;
		}
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 		
		Quaternion & Quaternion::operator*=( const Quaternion& _rQuat )
		{
		   REAL tmp[4];
	
		   tmp[0] = mfRealPart * _rQuat.mvImaginaryPart[0] + mvImaginaryPart[0] * _rQuat.mfRealPart + mvImaginaryPart[1] * _rQuat.mvImaginaryPart[2] - mvImaginaryPart[2] * _rQuat.mvImaginaryPart[1];
		   tmp[1] = mfRealPart * _rQuat.mvImaginaryPart[1] + mvImaginaryPart[1] * _rQuat.mfRealPart + mvImaginaryPart[2] * _rQuat.mvImaginaryPart[0] - mvImaginaryPart[0] * _rQuat.mvImaginaryPart[2];
		   tmp[2] = mfRealPart * _rQuat.mvImaginaryPart[2] + mvImaginaryPart[2] * _rQuat.mfRealPart + mvImaginaryPart[0] * _rQuat.mvImaginaryPart[1] - mvImaginaryPart[1] * _rQuat.mvImaginaryPart[0];
		   tmp[3] = mfRealPart * _rQuat.mfRealPart - mvImaginaryPart[0] * _rQuat.mvImaginaryPart[0] - mvImaginaryPart[1] * _rQuat.mvImaginaryPart[1] - mvImaginaryPart[2] * _rQuat.mvImaginaryPart[2];
	
		   mvImaginaryPart[0] = tmp[0];
		   mvImaginaryPart[1] = tmp[1];
		   mvImaginaryPart[2] = tmp[2];
		   mfRealPart = tmp[3];
		
		   return *this;
		} 

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 		
		Quaternion &  Quaternion::operator/=( REAL _fScalar ) 
		{
			mfRealPart /= _fScalar;
			mvImaginaryPart[0] /= _fScalar;
			mvImaginaryPart[1] /= _fScalar;
			mvImaginaryPart[2] /= _fScalar;
		
			return *this;
		}
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 		
		Quaternion &  Quaternion::operator*=(REAL _fScalar )	
		{
			mfRealPart *= _fScalar;
			mvImaginaryPart[0] *= _fScalar;
			mvImaginaryPart[1] *= _fScalar;
			mvImaginaryPart[2] *= _fScalar;
		
			return *this;
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		REAL Quaternion::GetLength( void ) const
		{
			return Math::SquareRoot	( mfRealPart * mfRealPart 
											+ mvImaginaryPart.GetX() * mvImaginaryPart.GetX() 
											+ mvImaginaryPart.GetY() * mvImaginaryPart.GetY() 
											+ mvImaginaryPart.GetZ() * mvImaginaryPart.GetZ() );
		}

		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		Quaternion& Quaternion::Normalize( void )
		{
			REAL fLen( GetLength() );

			mfRealPart		/= fLen;
			mvImaginaryPart /= fLen;

			return *this;
		}

		//---------------------------------------------------------------------
		// Dot product 
		//---------------------------------------------------------------------	
		REAL Quaternion::DotP( const Quaternion& _rQuatA, const Quaternion& _rQuatB )
		{
			return _rQuatA.ma[0] * _rQuatB.ma[0] + _rQuatA.ma[1] * _rQuatB.ma[1] + _rQuatA.ma[2] * _rQuatB.ma[2] + _rQuatA.ma[3] * _rQuatB.ma[3] ;
		}
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		void Quaternion::SLERP( REAL _fTime , const Quaternion & _rStart, const Quaternion & _rEnd, Quaternion & _rOutResult ) 
		{ 
			// ------------------------------
			// precondition 
			MB_ASSERT( Math::Abs( REAL(1.0f) - _rStart.GetLength() ) < REAL( 0.001f ) , "Should be unit length\n" );
			MB_ASSERT( Math::Abs( REAL(1.0f) - _rEnd.GetLength() ) < REAL( 0.001f ) , "Should be unit length\n" );
			// ------------------------------

			// Compute the cosine of the angle between the two vectors.
			REAL fDot = Quaternion::DotP(_rStart, _rEnd );

			Quaternion end;

			if( 0 <= fDot )
			{
				end = _rEnd;
			}
			else
			{
				end = -_rEnd;
				fDot = - fDot;
			}

			const double DOT_THRESHOLD = 0.9995f;
			if ( Math::Abs( fDot ) < DOT_THRESHOLD ) 
			{
				Math::Clamp( fDot, -1, 1 );				// Robustness: Stay within domain of acos()
				REAL fAngle = Math::ACos(fDot);  // theta_0 = angle between input vectors
  
				_rOutResult =	_rStart * Math::Sin( fAngle * ( 1 - _fTime ));
				_rOutResult +=	end * Math::Sin( fAngle * _fTime );
				_rOutResult /= Math::Sin( fAngle ); 
			}
			else
			{
				// If the inputs are too close for comfort, linearly interpolate
				// and normalize the result.
				NLERP(_fTime, _rStart, _rEnd, _rOutResult);
			}

			//return v0*cos(theta) + v2*sin(theta);

		}
		
		//---------------------------------------------------------------------------------------------- 
		//---------------------------------------------------------------------------------------------- 
		void Quaternion::NLERP( REAL _fTime , const Quaternion & _rStart, const Quaternion & _rEnd, Quaternion & _rOutResult ) 
		{
			// linearly interpolate each component, then normalize the Quaternion
			// Unlike spherical interpolation, this does not rotate at a constant velocity,
			// although that's not necessarily a bad thing
 
			_rOutResult = _rStart + (_rEnd - _rStart) * _fTime ;

			_rOutResult.Normalize();
		}		

	}
}



//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//	Implementations :	SRT 
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
namespace Melting 
{
	namespace Math
	{ 		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		SRT::SRT( void )
		{
			// /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !
		}
		
		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		SRT::SRT( const Vec3 & _rScale, const Vec3 & _rEulerAngles, const Vec3 & _rTranslate )
		:	mvScale			( _rScale )
		,	mquatRotation	( _rEulerAngles[0], _rEulerAngles[1], _rEulerAngles[2] )
		,	mvTranslate		( _rTranslate )
		{ 
		} 
 
		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		SRT::SRT( const Vec3 & _rScale, const Quaternion & _rquatRotation, const Vec3 & _rTranslate )
		:	mvScale			( _rScale )
		,	mquatRotation	( _rquatRotation )
		,	mvTranslate		( _rTranslate )
		{ 
		} 

		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		SRT::SRT( const SRT & _rSRT )
		:	mvScale			( _rSRT.mvScale ) 
		,	mquatRotation	( _rSRT.mquatRotation )	
		,	mvTranslate		( _rSRT.mvTranslate )
		{ 

		}			
		
		
		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		bool SRT::operator==( const SRT & _rOther )
		{
			return	mvScale == _rOther.mvScale 
				&&	mquatRotation == _rOther.mquatRotation 
				&&	mvTranslate == _rOther.mvTranslate;
		}
		
		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		bool SRT::operator!=( const SRT & _rOther )
		{
			return !( (*this) == _rOther );
		}
		

		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		void SRT::BuildIdentity( SRT& _rOutIdentity )
		{
			_rOutIdentity.mvScale = Vec3::ONE;
			_rOutIdentity.mquatRotation = Quaternion(1,0,0,0);
			_rOutIdentity.mvTranslate = Vec3::ZERO;
		}

		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		void SRT::BuildNull( SRT& _rOutNull )
		{
			_rOutNull.mvScale = Vec3::ZERO;
			_rOutNull.mquatRotation = Quaternion(0,0,0,0);
			_rOutNull.mvTranslate = Vec3::ZERO;
		}

		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		void SRT::Transform( const SRT& _rLocal, const SRT& _rParent, SRT& _rOutResult )
		{
			// TODO : OPTIMIZE 
			// scale
			_rOutResult.mvScale	= _rParent.mvScale * _rLocal.mvScale;

			// orientation
			_rOutResult.mquatRotation = _rParent.mquatRotation * _rLocal.mquatRotation ;

			// position 	 

			_rOutResult.mvTranslate =	_rParent.mquatRotation * ( _rLocal.mvTranslate * _rParent.mvScale); 
			_rOutResult.mvTranslate += 	_rParent.mvTranslate;  
		}
		
		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		void SRT::Multiply( REAL _fFactor, const SRT& _rBase, SRT& _rOutResult )
		{
			_rOutResult.mvTranslate		= _rBase.mvTranslate	* _fFactor;
			_rOutResult.mquatRotation	= _rBase.mquatRotation	* _fFactor;
			_rOutResult.mvScale			= _rBase.mvScale		* _fFactor;
		}
		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		void SRT::MultiplyInPlace( REAL _fFactor, SRT& _rOutResult )
		{
			_rOutResult.mvTranslate		*= _fFactor;
			_rOutResult.mquatRotation	*= _fFactor;
			_rOutResult.mvScale			*= _fFactor;
		}
		
		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		void SRT::Add( const SRT& _rA, const SRT& _rB, SRT& _rOutResult )
		{
			_rOutResult.mvTranslate		= _rA.mvTranslate	+ _rB.mvTranslate;
			_rOutResult.mquatRotation	= _rA.mquatRotation	+ _rB.mquatRotation;
			_rOutResult.mvScale			= _rA.mvScale		+ _rB.mvScale;
		}
		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		void SRT::AddInPlace( const SRT& _rA, SRT& _rOutResult )
		{
			_rOutResult.mvTranslate		+= _rA.mvTranslate;
			_rOutResult.mquatRotation	+= _rA.mquatRotation;
			_rOutResult.mvScale			+= _rA.mvScale;
		}
		
		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		void SRT::Interpolate( REAL _fTime, const SRT& _rStart, const SRT& _rEnd, SRT& _rOutResult )
		{
			// TODO : OPTIMIZE
 
			Vec3::LERP( _fTime, _rStart.mvScale, _rEnd.mvScale, _rOutResult.mvScale );

			Quaternion::SLERP( _fTime, _rStart.mquatRotation, _rEnd.mquatRotation, _rOutResult.mquatRotation );

			Vec3::LERP( _fTime, _rStart.mvTranslate , _rEnd.mvTranslate , _rOutResult.mvTranslate  );
		}
		
		//------------------------------------------------------------------------------ 
		//------------------------------------------------------------------------------
		void SRT::Inverse( const SRT& _rBase, SRT& _rOutResult )
		{
			// TODO : OPTIMIZE

			_rOutResult.mvTranslate[0]			= - _rBase.mvTranslate[0]; 
			_rOutResult.mvTranslate[1]			= - _rBase.mvTranslate[1]; 
			_rOutResult.mvTranslate[2]			= - _rBase.mvTranslate[2]; 
			Quaternion::Inverse( _rBase.mquatRotation, _rOutResult.mquatRotation );
			_rOutResult.mvScale[0]				= REAL(1.f)/ _rBase.mvScale[0]; 
			_rOutResult.mvScale[1]				= REAL(1.f)/ _rBase.mvScale[1]; 
			_rOutResult.mvScale[2]				= REAL(1.f)/ _rBase.mvScale[2]; 
		}					
			
	} // namespace Math
		
} // namespace Melting
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''


//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//	Implementations :	AABBOX - Axis-Aligned Bounding Box, represented by its Min and Max vertex (i.e. Min = 'LEFT-BOTTOM-FRONT', Max = 'RIGHT-TOP-REAR')	
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
namespace Melting
{
	namespace Math
	{ 			 	 
		//------------------------------------------------------------------------------
		// Default CTOR : /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !
		//------------------------------------------------------------------------------
		AABBox::AABBox( void )
		{}
		
		//------------------------------------------------------------------------------
		// CTOR
		//------------------------------------------------------------------------------
		AABBox::AABBox( const Vec3 & _rvMin, const Vec3 & _rvMax )
		: 	mvMin 	( _rvMin )
		,	mvMax 	( _rvMax )
		{	
		}
		
		//------------------------------------------------------------------------------
		// Copy-CTOR
		//------------------------------------------------------------------------------
		AABBox::AABBox( const AABBox & _rBox )
		: 	mvMin 	( _rBox.mvMin )
		,	mvMax 	( _rBox.mvMax )
		{	
		}
					
		//------------------------------------------------------------------------------
		// Setup the AABBox with the given 'min' position (bottom left front) of the box, and its volume (with absolute components required!)
		//------------------------------------------------------------------------------
		void AABBox::SetupFromMinPositionAndVolume( const Vec3 & _rvMinPosition, const Vec3 & _rvVolume )
		{
			//---------------------------------------
			// pre condition		
			MB_ASSERT( REAL( 0 ) < _rvVolume[0], "The X component of the volume must be strictly positive." );
			MB_ASSERT( REAL( 0 ) < _rvVolume[1], "The Y component of the volume must be strictly positive." );
			MB_ASSERT( REAL( 0 ) < _rvVolume[2], "The Z component of the volume must be strictly positive." );		
			//---------------------------------------
			
			mvMin = _rvMinPosition;
			mvMax = _rvMinPosition + _rvVolume;
		}
		
		//------------------------------------------------------------------------------
		// Setup the AABBox with the given 'center' position of the box, and its volume / 2 (with absolute components required!)
		//------------------------------------------------------------------------------
		void AABBox::SetupFromCenterPositionAndHalfVolume( const Vec3 & _rvCenterPosition, const Vec3 & _rvHalfVolume )
		{
			//---------------------------------------
			// pre condition		
			MB_ASSERT( REAL( 0 ) < _rvHalfVolume[0], "The X component of the volume must be strictly positive." );
			MB_ASSERT( REAL( 0 ) < _rvHalfVolume[1], "The Y component of the volume must be strictly positive." );
			MB_ASSERT( REAL( 0 ) < _rvHalfVolume[2], "The Z component of the volume must be strictly positive." );
			//---------------------------------------
		
			mvMin = _rvCenterPosition - _rvHalfVolume ;
			mvMax = _rvCenterPosition + _rvHalfVolume ;
		}
		
		//------------------------------------------------------------------------------
		// Get the Min...
		//------------------------------------------------------------------------------
		const Vec3&	AABBox::GetMin( void ) const
		{
			return mvMin;
		}
		
		//------------------------------------------------------------------------------
		// Get the Max !!!
		//------------------------------------------------------------------------------
		const Vec3&	AABBox::GetMax( void ) const
		{
			return mvMax;
		}

		
		//------------------------------------------------------------------------------
		// Set the Min...
		//------------------------------------------------------------------------------
		void AABBox::SetMin( const Vec3& _rvMin )
		{
			mvMin = _rvMin;
		}

		//------------------------------------------------------------------------------
		// Set the Max !!!
		//------------------------------------------------------------------------------
		void AABBox::SetMax( const Vec3& _rvMax ) 
		{
			mvMax = _rvMax;
		}

		//------------------------------------------------------------------------------
		// Get the the Center of the AABBox
		//------------------------------------------------------------------------------		
		void AABBox::GetCenter( Vec3& _rvResult ) const
		{
			_rvResult = (mvMax + mvMin);
			_rvResult *= REAL(0.5f);
		}
		
		//------------------------------------------------------------------------------
		// Get Half the height (Y) of this AABBox
		//------------------------------------------------------------------------------
		REAL AABBox::GetHalfHeight( void ) const
		{
			return ( mvMax[1] - mvMin[1] ) * REAL( 0.5f );			
		}

		//------------------------------------------------------------------------------
		// Get the 'Positive vertex' according to the given normal
		//------------------------------------------------------------------------------
		void AABBox::GetPositiveVertex( const Vec3& _rvNormal, Vec3& _rvResult ) const
		{
			_rvResult[0] = _rvNormal[0] >= REAL( 0 ) ? mvMax[0] : mvMin[0];
			_rvResult[1] = _rvNormal[1] >= REAL( 0 ) ? mvMax[1] : mvMin[1];
			_rvResult[2] = _rvNormal[2] >= REAL( 0 ) ? mvMax[2] : mvMin[2]; 
		}
		
		//------------------------------------------------------------------------------
		// Get the 'Negative vertex' according to the given normal
		//------------------------------------------------------------------------------
		void AABBox::GetNegativeVertex( const Vec3& _rvNormal, Vec3& _rvResult ) const
		{
			_rvResult[0] = _rvNormal[0] < REAL( 0 ) ? mvMax[0] : mvMin[0];
			_rvResult[1] = _rvNormal[1] < REAL( 0 ) ? mvMax[1] : mvMin[1];
			_rvResult[2] = _rvNormal[2] < REAL( 0 ) ? mvMax[2] : mvMin[2]; 
		}
		
		//------------------------------------------------------------------------------
		// Snap the given Position to the closest plane of the AABBox
		//------------------------------------------------------------------------------
		void AABBox::SnapToClosestPlane( bool _abSkippedDimensions[3], Vec3 & _rvOutSnaped ) const
		{
			//-----------------------------
			// precondition
			MB_ASSERT( !(_abSkippedDimensions[0] && _abSkippedDimensions[1] && _abSkippedDimensions[2] ) , "Skipping all dimensions is forbidden !\n" ) ;
			//-----------------------------
		
			UINT32 ulClosestMinMax 	= 0;
			UINT32 ulClosestDim 		= 0;
			
			REAL fClosestDistance 	(1000);
			REAL fTmpDistance;
			
			/// ATTENTION : trick !!!!!!!
			const Vec3* avMinMax = & mvMin; 
			
			for( UINT32 ulMinMax = 0; ulMinMax < 2; ulMinMax ++ )
			{
				for( UINT32 ulDim = 0; ulDim < 3; ulDim++ )
				{
					if( ! _abSkippedDimensions[ ulDim ] )
					{
						fTmpDistance = Math::Abs( avMinMax[ ulMinMax ][ ulDim ] - _rvOutSnaped[ ulDim ] );
						
						if( fTmpDistance < fClosestDistance )
						{
							fClosestDistance = fTmpDistance;
							
							ulClosestMinMax	= ulMinMax;
							ulClosestDim	= ulDim;
						}						
					}
				}
			}
			
			_rvOutSnaped[ ulClosestDim ] = avMinMax[ ulClosestMinMax ][ ulClosestDim ];
		}
		
		//------------------------------------------------------------------------------
		// Get the center of the model, and store it inside _rvOutCenter
		//------------------------------------------------------------------------------
		void AABBox::ComputeCenter( Vec3&  _rvOutCenter ) const
		{
			_rvOutCenter = 	mvMax + mvMin ;
			_rvOutCenter *= REAL(0.5f) ;
		};
		
		//------------------------------------------------------------------------------
		// Is the given point inside the AABBox
		//------------------------------------------------------------------------------
		bool AABBox::IsInside( const Vec3& _rvPoint ) const
		{
			return 	mvMin[0] <= _rvPoint[0] && _rvPoint[0] <= mvMax[0]
				&&	mvMin[1] <= _rvPoint[1] && _rvPoint[1] <= mvMax[1]
				&&	mvMin[2] <= _rvPoint[2] && _rvPoint[2] <= mvMax[2];
		}
		
		//------------------------------------------------------------------------------
		// operator : ==
		//------------------------------------------------------------------------------
		bool AABBox::operator == ( const AABBox & _rBox ) const
		{
			return mvMin == _rBox.mvMin && mvMax == _rBox.mvMax;
		}
		
		//------------------------------------------------------------------------------
		// operator : !=
		//------------------------------------------------------------------------------
		bool AABBox::operator != ( const AABBox & _rBox ) const
		{
			return ! ((*this) == _rBox);
		}
		
		//------------------------------------------------------------------------------
		// Compute the union of two AABBox, and store the result in the given reference parameter 'rOutBoxResult'
		//------------------------------------------------------------------------------
		void AABBox::Union( const AABBox & _rBoxA, const AABBox & _rBoxB, AABBox & _rOutBoxResult )
		{
			// might be faster with binary | operator ? kamo ne...			
			Vec3::GetVectorMin( _rBoxA.mvMin, _rBoxB.mvMin, _rOutBoxResult.mvMin );				
			Vec3::GetVectorMax( _rBoxA.mvMax, _rBoxB.mvMax, _rOutBoxResult.mvMax );
		}

		//------------------------------------------------------------------------------
		// Compute the union of a given AABBox and a point, and store the result in the given reference parameter 'rOutBoxResult'
		//------------------------------------------------------------------------------
		void AABBox::Union(const AABBox & _rBoxA, const Vec3 & _vPoint, AABBox & _rOutBoxResult)
		{			
			// might be faster with binary | operator ? kamo ne...
			Vec3::GetVectorMin( _rBoxA.mvMin, _vPoint, _rOutBoxResult.mvMin );				
			Vec3::GetVectorMax( _rBoxA.mvMax, _vPoint, _rOutBoxResult.mvMax ); 
		}
		 
	}
}
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''


//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//	Implementations :	BOUNDING SPHERE
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
namespace Melting
{
	namespace Math
	{ 			 
		//--------------------------------------------------------------------------------------------------------------------------
		// Default CTOR : /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !
		//--------------------------------------------------------------------------------------------------------------------------			
		BoundingSphere::BoundingSphere( void )		
		{}
		
		//--------------------------------------------------------------------------------------------------------------------------
		// CTOR
		//--------------------------------------------------------------------------------------------------------------------------
		BoundingSphere::BoundingSphere( const Vec3& _rvCenter, REAL _fRadius ) 
		:	mvCenter 			( _rvCenter )
		,	mfRadius			( _fRadius ) 
		{  }
		 
		//--------------------------------------------------------------------------------------------------------------------------
		// Copy-CTOR
		//--------------------------------------------------------------------------------------------------------------------------
		BoundingSphere::BoundingSphere( const BoundingSphere& _rSphere ) 
		:	mvCenter 	( _rSphere.mvCenter )
		,	mfRadius	( _rSphere.mfRadius ) 
		{
			
		}

		//--------------------------------------------------------------------------------------------------------------------------
		// Setup this BoundingSphere as the tightest sphere containing the given box
		//--------------------------------------------------------------------------------------------------------------------------
		void BoundingSphere::SetupFromAABBoxInside( const AABBox & rBox )
		{  
			rBox.ComputeCenter( mvCenter );	// the center will change when the clip matrix will be applied (that is if the sphere is in the local coords)
			
			Vec3 vfxDiff = rBox.GetMax() - rBox.GetMin();
			mfRadius		= vfxDiff.GetLength() / REAL(2.f);
		}
		
		//--------------------------------------------------------------------------------------------------------------------------
		// Setup this BoundingSphere as the biggest sphere contained in the given box
		//--------------------------------------------------------------------------------------------------------------------------
		void BoundingSphere::SetupFromAABBoxOutside( const AABBox & rBox )
		{ 
			rBox.ComputeCenter( mvCenter );	// the center will change when the clip matrix will be applied (that is if the sphere is in the local coords)
			
			Vec3 size 		= rBox.GetMax() - rBox.GetMin(); 
			mfRadius		= min( min( size[0], size[1] ), size[2] ) * REAL(0.5f);	 
		}

		//--------------------------------------------------------------------------------------------------------------------------
		// Get the BoundingSphere center
		//--------------------------------------------------------------------------------------------------------------------------
		const Vec3& BoundingSphere::GetCenter( void ) const
		{
			return mvCenter;
		}
		
		//--------------------------------------------------------------------------------------------------------------------------
		// Set the BoundingSphere center
		//--------------------------------------------------------------------------------------------------------------------------
		void BoundingSphere::SetCenter( const Vec3& _rvCenter )
		{
			mvCenter = _rvCenter;
		}
		 
		//--------------------------------------------------------------------------------------------------------------------------
		// Get the BoundingSphere radius
		//--------------------------------------------------------------------------------------------------------------------------
		REAL BoundingSphere::GetRadius( void ) const
		{
			return mfRadius;
		}
		
		//--------------------------------------------------------------------------------------------------------------------------
		// Set the BoundingSphere radius
		//--------------------------------------------------------------------------------------------------------------------------
		void  BoundingSphere::SetRadius( REAL _fRadius ) 
		{
			mfRadius = _fRadius; 
		}		
			

		//--------------------------------------------------------------------------------------------------------------------------
		// Is the given point inside the BoundingSphere
		//--------------------------------------------------------------------------------------------------------------------------
		bool BoundingSphere::IsInside( const Vec3& _rvPoint ) const
		{
			return ( _rvPoint - mvCenter ).GetSquaredLength() < ( mfRadius * mfRadius );
		}
			
		//--------------------------------------------------------------------------------------------------------------------------
		// operator : ==
		//--------------------------------------------------------------------------------------------------------------------------
		bool BoundingSphere::operator == ( const BoundingSphere & _rBoundingSphere ) const	
		{ 
			return mfRadius == _rBoundingSphere.mfRadius && mvCenter == _rBoundingSphere.mvCenter;
		}
		
		//--------------------------------------------------------------------------------------------------------------------------
		// operator : !=
		//--------------------------------------------------------------------------------------------------------------------------
		bool BoundingSphere::operator != ( const BoundingSphere & _rBoundingSphere ) const
		{ 
			return !(( *this ) == _rBoundingSphere);	
		} 
	 
//#include <nitro/itcm_begin.h>
		//--------------------------------------------------------------------------------------------------------------------------
		// Compute the union of rSphereA and rSphereB, and store the result in rOutUnionSphere -- PS : rSphereA or rSphereB can be used as rOutUnionSphere 
		//--------------------------------------------------------------------------------------------------------------------------
		void BoundingSphere::Union( const BoundingSphere & _rSphereA , const BoundingSphere & _rSphereB, BoundingSphere & _rOutUnionSphere )
		{ 
			//*/
			// vector from a center to the other
			Vec3 vCenterToCenter 	= (_rSphereB.mvCenter - _rSphereA.mvCenter); 
			REAL fSqrDistance 		= vCenterToCenter.GetSquaredLength();			
			REAL fSqrRadiusDiff 	= Math::Square( _rSphereB.GetRadius() - _rSphereA.GetRadius() );
			
	        // The sphere with the larger radius encloses the other;
	        // just set s to be the larger of the two spheres
			if( fSqrRadiusDiff >= fSqrDistance )
			{
				if( _rSphereB.GetRadius() >= _rSphereA.GetRadius() )
				{
					_rOutUnionSphere = _rSphereB;
				}
				else
				{
					_rOutUnionSphere = _rSphereA;
				}
			}
	        // Spheres partially overlapping or disjoint
			else
			{
				REAL fDistance = Math::SquareRoot( fSqrDistance );
				_rOutUnionSphere.SetRadius( (fDistance + _rSphereB.GetRadius() + _rSphereA.GetRadius()) * REAL(0.5f) ) ;
				
				_rOutUnionSphere.SetCenter( _rSphereA.GetCenter() );
				
				const REAL EPSILON ( 0.01f );
				
				if( fDistance > EPSILON )
				{
					Vec3 vNewCenterOffset = vCenterToCenter * ( ( _rOutUnionSphere.GetRadius() - _rSphereA.GetRadius() ) / fDistance );
					_rOutUnionSphere.SetCenter( _rOutUnionSphere.GetCenter() + vNewCenterOffset );
				}				
			}
			/*/
			
			// vector from a center to the other
			Vec3 vCenterToCenter = (_rSphereA.mvCenter - _rSphereB.mvCenter); 
			
			// if the 2 centers are not too close
			if( Vec3( REAL( 0 )) != vCenterToCenter ) 
			{ 
				// distance² between the 2 centers
				REAL fxSqrDistance = vCenterToCenter.GetSquaredLength();
				
				// 
				REAL fxSqrRadiusDiff = (_rSphereA.GetRadius() - _rSphereB.GetRadius()).Squared();
				
				// if one is completely inside another
				if( fxSqrRadiusDiff >= fxSqrDistance )
				{
					// if B radius is greater than A's, then A is inside B
					if( _rSphereA.GetRadius() < _rSphereB.GetRadius())
					{
						_rOutUnionSphere = _rSphereB;
						return;
					}
					else
					{
						_rOutUnionSphere = _rSphereA;
						return;					
					}
				}
				// 
				else
				{	
					vCenterToCenter.Normalize();
					
					Vec3 vEdgeA 	= _rSphereA.mvCenter;
					vEdgeA 			+= vCenterToCenter * _rSphereA.mfRadius;
					
					Vec3 vEdgeB 	= _rSphereB.mvCenter;
					vEdgeB			-= vCenterToCenter * _rSphereB.mfRadius;
					
					// radius = Length ( A - B ) / 2
					_rOutUnionSphere.mfRadius 	= (vEdgeA - vEdgeB).GetLength() ;
					_rOutUnionSphere.mfRadius 	*= REAL(0.5f);	
					// center = ( A + B ) / 2
					_rOutUnionSphere.mvCenter	= (vEdgeA + vEdgeB) ;
					_rOutUnionSphere.mvCenter	*= REAL(0.5f);  
				}
			}
			else
			{
				// if B radius is greater than A's, then A is inside B
				if( _rSphereA.GetRadius() < _rSphereB.GetRadius())
				{
					_rOutUnionSphere = _rSphereB;
					return;
				}
				else
				{
					_rOutUnionSphere = _rSphereA;
					return;					
				}	
			} 
			//*/
		}
//#include <nitro/itcm_end.h>
	}
}

//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//	Implementations :	CONE 
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
namespace Melting
{
	namespace Math
	{ 			  
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		Cone::Cone( void )  
		{ 
			// /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !		
		}
		
		//------------------------------------------------------------------------------
		// CTOR
		//------------------------------------------------------------------------------
		Cone::Cone( const Vec3& _rvVertex, const Vec3& _rvDirection, REAL fLength, REAL _fAngle )
		:	mvVertex		( _rvVertex )
		,	mvDirection		( _rvDirection )
		,	mfLength		( fLength )
		, 	mfAngle		( _fAngle ) 
		{
			ComputeSinCos();
			mvDirection.Normalize();
		}
		
		//------------------------------------------------------------------------------
		// Copy CTOR
		//------------------------------------------------------------------------------
		Cone::Cone( const Cone& _rCone )
		:	mvVertex		( _rCone.mvVertex )
		,	mvDirection		( _rCone.mvDirection )
		,	mfLength		( _rCone.mfLength )
		, 	mfAngle			( _rCone.mfAngle )
		,	mfHalfAngleCos	( _rCone.mfHalfAngleCos )
		,	mfHalfAngleSin	( _rCone.mfHalfAngleSin )
		{
		}
		
		//------------------------------------------------------------------------------
		// Setup the cone from a starting point (the vertex), and endind point, and an angle
		//------------------------------------------------------------------------------
		void Cone::SetupFromTwoPointsAndAngle( const Vec3& _rvStart, const Vec3& _rvEnd, REAL _fAngle )
		{ 
			Vec3 vDirection = _rvEnd - _rvStart;
			
			SetLength	( vDirection.GetLength() );
			SetDirection( vDirection.Normalize() ); 
			SetVertex	( _rvStart );
			SetAngle	( _fAngle );
		}

		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		void Cone::SetVertex( const Vec3& _rvVertex )  	 		
		{ 
			mvVertex = _rvVertex ; 
		}
		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		const Vec3& Cone::GetVertex( void ) const 							
		{ 
			return mvVertex; 
		}
		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		void Cone::SetDirection( const Vec3& _rvDirection )  		 
		{ 		
			// ----------------------------
			// precondition
			ASSERT_UNIT_LENGTH_VECTOR( _rvDirection );
			// ----------------------------
				
			mvDirection = _rvDirection; 
		}		
		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		const Vec3& Cone::GetDirection( void ) const 								
		{ 
			return mvDirection; 
		}
		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		void Cone::SetLength( REAL _fLength ) 		 						
		{ 
			mfLength = _fLength;
		}
		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		REAL Cone::GetLength( void ) const	 		 						
		{ 
			return mfLength; 
		}
		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		void Cone::SetAngle( REAL _fAngle ) 		 						
		{ 
			mfAngle = _fAngle; 
			ComputeSinCos(); 
		}
		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		REAL Cone::GetAngle( void ) const	 		 						
		{ 
			return mfAngle; 
		}
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		REAL Cone::GetHalfAngleCos( void ) const							
		{ 
			return mfHalfAngleCos; 
		}
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		REAL Cone::GetHalfAngleSin( void ) const							
		{ 
			return mfHalfAngleSin; 
		} 		 		 
		
		//------------------------------------------------------------------------------
		// Is the Sphere inside the Cone ??
		//------------------------------------------------------------------------------
		bool Cone::IsInside( const BoundingSphere & _rSphere, const Cone & _rCone )
		{	 
			// ----------------------------
			// precondition
			ASSERT_UNIT_LENGTH_VECTOR( _rCone.GetDirection() );
			// ----------------------------
			
			// 1. we offset the cone vertex P to the rear from sphere radius units, so that testing if the sphere center O is inside the 
			//		'offseted cone' will be equivalent to 'the sphere is inside the original cone' 
			//		(...modulo a small error when the sphere center is actually behind the original cone vertex... but that's enough for our needs
			// 2. we project the sphere center O against the cone direction and get the length of this projected vector as PO'
			// 		-> if ( length (PO') > length of the vector ... ) then it's inside the cone !
							
			Cone cone ( _rCone.GetVertex() - _rCone.GetDirection() * _rSphere.GetRadius(), _rCone.GetDirection(), _rCone.GetLength(), _rCone.GetAngle() )  ;
			
			return Cone::IsInside( _rSphere.GetCenter() , cone ); 
		}	
		
		//------------------------------------------------------------------------------
		// Is the Point inside the Cone ??
		//------------------------------------------------------------------------------
		bool Cone::IsInside( const Vec3& _rPoint, const Cone& _rCone )
		{
			// ----------------------------
			// precondition
			ASSERT_UNIT_LENGTH_VECTOR(_rCone.GetDirection());
			// ----------------------------
										
			// we project the point P against the cone direction and get the length of this projected vector as PO'
			// 		-> if ( length (PO') > length of the vector ... ) then it's inside the cone !				 
			Vec3 vVertexToCenter = _rPoint - _rCone.GetVertex() ;
			
			if( Vec3( REAL( 0 )) != vVertexToCenter )
			{ 		
//				LOG( "Point :\n" ); _rPoint.Print();
//				LOG( "Cone vertex:\n" ); _rCone.GetVertex().Print();
//				LOG( "vertex to point:\n" ); vVertexToCenter.Print();
				
				REAL fDistVertexToProjectedCenter	= _rCone.GetDirection().DotP( vVertexToCenter );
				
				// if the projected center is not beyond the length, the point can be inside 
				if( _rCone.GetLength() >= fDistVertexToProjectedCenter )
				{
					REAL fDistMinToInsideCone	= vVertexToCenter.GetLength() ;
					fDistMinToInsideCone		*= _rCone.GetHalfAngleCos();
					
					return fDistVertexToProjectedCenter >= fDistMinToInsideCone;	
				}
				// if the projected center is beyond the length, the point IS NOT inside anyway
				else
				{
					return false;
				}
			}
			// the point IS at the same position than the cone vertex, so it IS inside
			else
			{
				return true;
			} 
		}
		
		//------------------------------------------------------------------------------
		// Compute the Sin & Cos of the (Angle/2) 
		//------------------------------------------------------------------------------	
		void Cone::ComputeSinCos( void )
		{
			const REAL fHal_fAngle = mfAngle * REAL(0.5f);
			
			mfHalfAngleCos = Math::Cos( fHal_fAngle );
			mfHalfAngleSin = Math::Sin( fHal_fAngle );
		}
		
		
		//------------------------------------------------------------------------------	
		// operator : == 
		//------------------------------------------------------------------------------	
		bool Cone::operator == ( const Cone & _rCone ) const
		{
			return 		mvVertex 	== _rCone.mvVertex 
					&& 	mvDirection == _rCone.mvDirection 
					&& 	mfLength 	== _rCone.mfLength 
					&& 	mfAngle 	== _rCone.mfAngle;	
		}
		
		//------------------------------------------------------------------------------	
		// operator : != 
		//------------------------------------------------------------------------------	
		bool Cone::operator != ( const Cone & _rCone ) const
		{
			return !(*this == _rCone);	
		}
	}
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//	Implementations :	PLANE 
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
namespace Melting
{
	namespace Math
	{ 
		//------------------------------------------------------------------------------		
		// Default CTOR : /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !
		//------------------------------------------------------------------------------		
		Plane::Plane ( void )
		{}
		
		//------------------------------------------------------------------------------		
		// CTOR
		//------------------------------------------------------------------------------		
		Plane::Plane ( const Vec3 & _rvNormal , REAL _fDistanceToOrigin  )
		:	mvNormal 	( _rvNormal )
		,	mfDist		( - _fDistanceToOrigin )
		{
			// ----------------------------
			// postcondition
			ASSERT_UNIT_LENGTH_VECTOR( mvNormal );
			// ----------------------------
		}
		
		//------------------------------------------------------------------------------		
		// Copy - CTOR
		//------------------------------------------------------------------------------		
		Plane::Plane ( const Plane & _rPlane )
		:	mvNormal 	( _rPlane.mvNormal )
		,	mfDist		( _rPlane.mfDist )
		{
			// ----------------------------
			// postcondition
			ASSERT_UNIT_LENGTH_VECTOR( mvNormal );
			// ----------------------------
		} 
		//------------------------------------------------------------------------------		
		// Create a Plane defined by a triangle
		//------------------------------------------------------------------------------			
		void Plane::SetupFromTriangle(const Vec3 & _rvVertex1, const Vec3 & _rvVertex2, const Vec3 & _rvVertex3)
		{ 
			Vec3 v1 = _rvVertex2 - _rvVertex1;
			Vec3 v2 = _rvVertex3 - _rvVertex1; 
 			
			SetupFromNormalAndPoint( v1.CrossP(v2).Normalize(), _rvVertex1); 
		}


		//------------------------------------------------------------------------------
		// Get the Plane normal
		//------------------------------------------------------------------------------
		const Vec3& Plane::GetNormal( void ) const
		{
			return mvNormal;
		}
		//------------------------------------------------------------------------------
		// Get the Plane distance to the world origin (along the direction)
		//------------------------------------------------------------------------------
		REAL Plane::GetDistance( void ) const
		{
			return mfDist;	
		}
		

//#include <nitro/itcm_begin.h>
		//------------------------------------------------------------------------------		
		// Create a plane defined by a normal and a point on this plane
		//------------------------------------------------------------------------------
		void Plane::SetupFromNormalAndPoint( const Vec3 & _rvNormal, const Vec3 & _rvVertex )
		{
			// ----------------------------
			// precondition
			ASSERT_UNIT_LENGTH_VECTOR( _rvNormal );
			// ----------------------------
		
			mvNormal 	= _rvNormal;
			mfDist 		= - mvNormal.DotP( _rvVertex );
		} 
		
		//------------------------------------------------------------------------------
		// Compute the distance from the point to the plane 
		//------------------------------------------------------------------------------
		REAL Plane::DistanceFrom( const Vec3 & _rvVertex ) const
		{
			return ( mvNormal.DotP( _rvVertex ) + mfDist ) ;
		} 
		
		//------------------------------------------------------------------------------
		// Project the given point 'rvVertex 'on this Plane, and store the _result in 'rvOutProjectedVertex';
		//------------------------------------------------------------------------------
		void Plane::Project( const Vec3 & _rvVertex , Vec3 & _rvOutProjectedVertex )
		{
			_rvOutProjectedVertex = _rvVertex - mvNormal * DistanceFrom( _rvVertex );
		} 

		//------------------------------------------------------------------------------
		// Tells wether the given point is above this Plane -- /!\ Assume the normal points toward the "ABOVE"  or "FRONT" 
		//------------------------------------------------------------------------------
		bool Plane::IsAbove(const Vec3 & _rvVertex ) const
		{ 
			return ( DistanceFrom( _rvVertex ) > REAL( 0 ) ) ;
		} 
	
//#include <nitro/itcm_end.h>	
		//------------------------------------------------------------------------------
		// operator : == 
		//------------------------------------------------------------------------------			
		bool Plane::operator == ( const Plane & _rPlane ) const
		{
			return mvNormal == _rPlane.mvNormal && mfDist == _rPlane.mfDist;
		}
		
		//------------------------------------------------------------------------------
		// operator : != 
		//------------------------------------------------------------------------------
		bool Plane::operator != ( const Plane & _rPlane ) const
		{
			return !((*this) == _rPlane );
		}
		
		
	} // namespace Math
		
} // namespace Melting
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
 
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
//	Implementations :	RAY 
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
namespace Melting
{
	namespace Math
	{ 
		//------------------------------------------------------------------------------
		// Default CTOR : 
		//------------------------------------------------------------------------------
		Ray::Ray( void )
		{
			// WARNING : /!\ Keep it uninitialized for performance - cause using a default CTOR means we will initialize it differently later !
		}
		
		//------------------------------------------------------------------------------
		// CTOR - the Direction must be unit length
		//------------------------------------------------------------------------------
		Ray::Ray( const Vec3 & _rvOrigin, const Vec3 & _rvDirection )
		:	mvOrigin			( _rvOrigin )
		,	mvDirection			()
		,	mvInverseDirection	()
		,	mvSignedNeg			()
		{
			SetDirection( _rvDirection );
		
			// ----------------------------
			// postcondition
			// NDAL : already tested within    SetDirection()
			// ASSERT_UNIT_LENGTH_VECTOR(mvDirection);
			// ----------------------------
		} 
	
		//------------------------------------------------------------------------------
		// Copy CTOR
		//------------------------------------------------------------------------------
		Ray::Ray( const Ray & _rRay )
		:	mvOrigin			( _rRay.mvOrigin )
		,	mvDirection			( _rRay.mvDirection )
		,	mvInverseDirection	( _rRay.mvInverseDirection )
		,	mvSignedNeg			( _rRay.mvSignedNeg )
		{
			// ----------------------------
			// postcondition
			ASSERT_UNIT_LENGTH_VECTOR(mvDirection);
			// ----------------------------
		}
	
		//------------------------------------------------------------------------------	
		// Create a Ray defined by a segment composed of 2 points : Start & End
		//------------------------------------------------------------------------------
		void Ray::SetupFromSegment(const Vec3 & _rvVertexStart, const Vec3 & _rvVertexEnd )
		{
			SetOrigin	( _rvVertexStart );
			// FIXME !!!!!
			SetDirection( ( _rvVertexEnd - _rvVertexStart ).Normalize() ); 
		}
		
		
		//------------------------------------------------------------------------------	
		// Set the origin of the Ray
		//------------------------------------------------------------------------------	
		void Ray::SetOrigin( const Vec3 & _rvOrigin ) 
		{
			mvOrigin = _rvOrigin;
		}
		//------------------------------------------------------------------------------	
		// Get the origin of the Ray
		//------------------------------------------------------------------------------	
		const Vec3 & Ray::GetOrigin( void ) const
		{
			return mvOrigin;
		} 

		//------------------------------------------------------------------------------	
		// Set the direction of the Ray
		//------------------------------------------------------------------------------	
		void Ray::SetDirection( const Vec3 & _rvDirection )
		{
			// ----------------------------
			// precondition
			ASSERT_UNIT_LENGTH_VECTOR( _rvDirection );
			// ----------------------------
		
			mvDirection = _rvDirection;
			Vec3::GetVectorInverse( mvDirection, mvInverseDirection );
			
			mvSignedNeg[ 0 ] = (mvDirection[ 0 ] < REAL(0)) ? 1 : 0;
			mvSignedNeg[ 1 ] = (mvDirection[ 1 ] < REAL(0)) ? 1 : 0;
			mvSignedNeg[ 2 ] = (mvDirection[ 2 ] < REAL(0)) ? 1 : 0;
		}
			
		//------------------------------------------------------------------------------	
		// Get the direction of the Ray
		//------------------------------------------------------------------------------	
		const Vec3 & Ray::GetDirection( void ) const
		{
			return mvDirection;			
		}
		//------------------------------------------------------------------------------	
		// Get the 1/direction of the Ray
		//------------------------------------------------------------------------------	
		const Vec3 & Ray::GetInverseDirection( void ) const
		{
			return mvInverseDirection;			
		}

//#include <nitro/itcm_begin.h>
		//------------------------------------------------------------------------------	
		// Get the 3D position of a point on the Ray, located at the given distance from the Ray origin -- the _result is stored in the _reference parameter 'rvOutPoint'
		//------------------------------------------------------------------------------	
		void Ray::ComputePoint( REAL _fDistanceFromOrigin, Vec3 & _rvOutPoint ) const
		{
			_rvOutPoint = mvOrigin + mvDirection * _fDistanceFromOrigin;
		}

		//------------------------------------------------------------------------------
		// Compute the distance from the Ray origin to the given Plane
		//------------------------------------------------------------------------------
		REAL Ray::ComputeDistance( const Plane & _rPlane ) const
		{			
			// ----------------------------
			// precondition
			ASSERT_UNIT_LENGTH_VECTOR( _rPlane.GetNormal() ) ;
			// ----------------------------
		
			REAL fxCosAngle = _rPlane.GetNormal().DotP( mvDirection );
			
			MB_ASSERT( REAL(0) != fxCosAngle, "Ray is parallel to the plane, so its distance is infinity!" );
			
			REAL fxInv = Math::Inverse( fxCosAngle );
			
			return - _rPlane.DistanceFrom( mvOrigin ) * fxInv;				
		}
			
		//------------------------------------------------------------------------------
		// Compute the intersection of the given Plane with this Ray , and store the _result in 'rvOutIntersect', and _returns wether there was an intersection or not
		//------------------------------------------------------------------------------
		bool Ray::ComputeIntersection( const Plane & _rPlane, Vec3 & _rvOutIntersect ) const
		{ 			
			// ----------------------------
			// precondition
			ASSERT_UNIT_LENGTH_VECTOR( _rPlane.GetNormal() );
			// ----------------------------
		
			// Ray-Plane intersection !
			// vRay ( vOrig, vDir) : 		vPosition = vOrig + vDir * T
			// Plane ( vNormal, fxDist) : 	vPosition <dot> vNormal + fxDist;
			// obviously, we need T, 
			// substituting for vPosition, we have 
			// T = - ( vOrig <dot> vNormal + fxDist) / ( vNormal <dot> vDir )
			// where vNormal < dot > vDir is also the Cosine of the angle between these 2 vectors
			
			bool bIntersect = false;
				
			REAL fxCosAngle = _rPlane.GetNormal().DotP( mvDirection );
			
			if( REAL( 0 ) != fxCosAngle )
			{ 
				REAL fxDistanceAlongRay = - _rPlane.DistanceFrom( mvOrigin ) / fxCosAngle ;
				
				ComputePoint( fxDistanceAlongRay, _rvOutIntersect );
				
				bIntersect = true;
			}
			
			return bIntersect;
		}
 
		//------------------------------------------------------------------------------
		// Returns wether this Ray intersects the given Plane or not
		//------------------------------------------------------------------------------
		bool Ray::DoesIntersect( const Plane & _rPlane ) const
		{ 			
			// ----------------------------
			// precondition
			ASSERT_UNIT_LENGTH_VECTOR( _rPlane.GetNormal() );
			// ----------------------------
			
			return REAL( 0 ) != _rPlane.GetNormal().DotP( GetDirection() );
		} 

		//--------------------------------------------------------------------------------------------------------------------------
		// Returns wether this Ray intersects the given AABBox or not
		//--------------------------------------------------------------------------------------------------------------------------
		bool Ray::DoesIntersect( const AABBox & _rAABBox ) const
		{  
			Vec3 vIntersection ( REAL( 0 ) );
			REAL intersectionDistance ( 0 );
			return ( ComputeIntersection( _rAABBox , vIntersection , intersectionDistance ) ); 		
		}
		
		//--------------------------------------------------------------------------------------------------------------------------
		// Intersect _ray R(t) = p + t*d against AABB a. When intersecting,
		// _return intersection distance tmin and point q of intersection
		//--------------------------------------------------------------------------------------------------------------------------
		bool Ray::ComputeIntersection( const AABBox & _rAABBox, Vec3& _rvOutIntersection, REAL& _rOutIntersectionDistance ) const	
		{		
		    _rOutIntersectionDistance = REAL( 0 );	// set to -REAL_MAX to get first hit on line
		    REAL tmax ( 500000 ); 					// set to REAL_MAX max distance _ray can travel (for segment)			 
 		
		    // For all three slabs ( SLAB == space between 2 parallel planes )
		    for ( UINT32 ulDim = 0; ulDim < 3; ulDim++ )
		    {
		        if ( REAL(0) == mvDirection[ulDim] ) 
		        {
		            // Ray is parallel to slab. No hit if origin not within slab
		            if (	mvOrigin[ulDim] < _rAABBox.GetMin()[ulDim] 
		            	|| 	mvOrigin[ulDim] > _rAABBox.GetMax()[ulDim] )
		            {
		            	return false;	
		            }
		        } 
		        else 
		        {
		            // Compute intersection t value of _ray with near and far plane of slab
		            REAL t1 = ( _rAABBox.GetMin()[ ulDim ] - mvOrigin[ ulDim ] ) * mvInverseDirection[ ulDim ];
		            REAL t2 = ( _rAABBox.GetMax()[ ulDim ] - mvOrigin[ ulDim ] ) * mvInverseDirection[ ulDim ];
		            
		            // Make t1 be intersection with near plane, t2 with far plane
		            if (t1 > t2) 
		            {
		            	Math::Swap(t1, t2);	
		            }
		            // Compute the intersection of slab intersections intervals
		            _rOutIntersectionDistance = max(_rOutIntersectionDistance, t1);
		            
		            tmax = min(tmax, t2);
		            
		            // Exit with no collision as soon as slab intersection becomes empty
		            if (_rOutIntersectionDistance > tmax) 
		            {		       
		            	return false;	
		            }
		        }
		    }
		    
		    // Ray intersects all 3 slabs. Return point (q) and intersection t value (tmin) 
		    _rvOutIntersection = mvOrigin + mvDirection * _rOutIntersectionDistance;
		    
		    return true;
		}
		
		//--------------------------------------------------------------------------------------------------------------------------
		// same as ComputeIntersection, except that we can specify the ScaleRotateTranslate of the box !
		//--------------------------------------------------------------------------------------------------------------------------
		bool Ray::ComputeIntersection ( const AABBox & _rAABBox, const SRT& _rBoxSRT, Vec3& _rvOutIntersection, REAL& _rOutIntersectionDistance ) const
		{ 
			// -> on passe le _rayon dans le _repère de la box !!!!
			AABBox box ( _rAABBox.GetMax() * _rBoxSRT.mvScale , _rAABBox.GetMax() * _rBoxSRT.mvScale );

			Quaternion quatInverse = _rBoxSRT.mquatRotation.GetInverse();

			Vec3 vRayOrigin = this->GetOrigin();
			vRayOrigin -= _rBoxSRT.mvTranslate;			
			vRayOrigin = quatInverse.Rotate( vRayOrigin );

			Vec3 vDirection = quatInverse.Rotate( this->GetDirection() );
			Ray _ray( vRayOrigin, vDirection );

			return _ray.ComputeIntersection( box, _rvOutIntersection, _rOutIntersectionDistance );			
		}

		//--------------------------------------------------------------------------------------------------------------------------
		// Returns wether this Ray intersects the given BoundingSphere or not
		//--------------------------------------------------------------------------------------------------------------------------
		bool Ray::DoesIntersect( const BoundingSphere & _rBSphere ) const
		{ 
			Vec3 	vRayOriginToSphereCenter 				= _rBSphere.GetCenter() - mvOrigin;
			REAL	fDistanceFromOriginToProjectedCenter 	= mvDirection.DotP( vRayOriginToSphereCenter );
			
			// if the sphere center is completely in FRONT of the Origin 
			if	( 	REAL( 0 ) < fDistanceFromOriginToProjectedCenter 
				&& 	Math::Square( _rBSphere.GetRadius() ) < ( fDistanceFromOriginToProjectedCenter * fDistanceFromOriginToProjectedCenter ) 
				)
			{
				Vec3 	vProjectedCenter;
				ComputePoint( fDistanceFromOriginToProjectedCenter, vProjectedCenter );
				
				REAL 	fSqrDistanceFromProjectedCenterToCenter = (vProjectedCenter - _rBSphere.GetCenter()).GetSquaredLength();
				
				// if the distance from projectedCenter to the _ray is greater than the _ray then there is no intersection
				return ( fSqrDistanceFromProjectedCenterToCenter < Math::Square( _rBSphere.GetRadius() ) );
			}
			// if the sphere is BEHIND the origin (or not completely in front of it)
			else
			{
				return false;
			}
			
			/*
			// TODO : another version :: check if it's faster!!!
			
			// Assumes _ray.d is normalized
			// Uses quadratic formula (AX² + BX + C = 0), A = 1 because _ray.d is normalized.
			// For derivation, see:
			// http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter1.htm			
			Vec3 vDiff = GetOrigin() - _rBSphere.GetCenter();
			
			REAL fxB = GetDirection().DotP( vDiff ); 					// B _really equals 2.0f this value
			REAL fxC = vDiff.DotP( vDiff ) - ( _rBSphere.GetRadius() * _rBSphere.GetRadius()) ;	// Discriminant.  D _really equals 4.0f times this value
			REAL fxDiscriminant = fxB * fxB - fxC;							// the common ( B² - 4AC )
			
			return fxDiscriminant > REAL( 0 ); 
			*/
		} 	  
		
		//--------------------------------------------------------------------------------------------------------------------------
		// operator : ==
		//--------------------------------------------------------------------------------------------------------------------------
		bool Ray::operator == ( const Ray & _rRay ) const
		{
			return mvOrigin == _rRay.mvOrigin && mvDirection == _rRay.mvDirection;
		}
			
		//--------------------------------------------------------------------------------------------------------------------------
		// operator : !=
		//--------------------------------------------------------------------------------------------------------------------------
		bool Ray::operator != ( const Ray & _rRay ) const
		{
			return !((*this) == _rRay );
		}
					
			
	} // namespace Math
		
} // namespace Melting
 
//''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

//#include "UnitTest.h"

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//														UNIT TESTs
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Melting
{
	namespace Math
	{
		// ------------------------------------------------------------------------------
		// ------------------------------------------------------------------------------
		bool UtilsUnitTest( void )
		{			
			LOG("::::::::: Functions : Math from Utils.h ::::::::::\n");
	
			UTEST_START; 
			// --- -2. --------------------------------------------------------------------
			LOG("-- Math::Equal( REAL a, REAL b )\n");
			{ 
				UTEST( Math::Equal( 1.0f, 1.0f ) );		 
				UTEST( Math::Equal( 1.0f, 1.0f + 0.00001f, 0.0001f ) ); 
				UTEST( Math::Equal( 1.0f, 1.0f - 0.00001f, 0.0001f ) ); 
			}
			// --- -1. --------------------------------------------------------------------
			LOG("-- Math::Abs( REAL ) \n");
			{	
				UTEST( 0.0f		== Math::Abs( 0.0f ) );	
				UTEST( 0.0f		== Math::Abs( -0.0f ) );
				UTEST( 1.0f		== Math::Abs( 1.0f ) );	
				UTEST( 1.0f		== Math::Abs( - 1.0f ) );	
				UTEST( 12.0f	== Math::Abs( 12.0f ) );	
				UTEST( 12.0f	== Math::Abs( - 12.0f ) );	
			}
			// --- 0. --------------------------------------------------------------------
			LOG("-- Math::Inverse( REAL )\n");
			{ 
				UTEST( 1.0f 	== Math::Inverse( 1.0f ) );	
				UTEST( 2.0f 	== Math::Inverse( 0.5f ) );	
				UTEST( 0.5f 	== Math::Inverse( 2.0f ) );		 
			}
			// --- 1. --------------------------------------------------------------------
			LOG("-- Math::Swap( REAL, REAL )\n");
			{ 
				REAL A = 2.0f;
				REAL B = 0.0f;

				REAL expectedA = 0.0f;
				REAL expectedB = 2.0f;

				Math::Swap( A, B );

				UTEST( expectedA == A && expectedB == B );	 
			}
			// --- 2. --------------------------------------------------------------------
			LOG("-- Math::Clamp( REAL & _rfOutClamped, REAL _fMin, REAL _fMax )\n");
			{
				{	
					REAL fValue = 10.f;
					Math::Clamp( fValue, 0.0f, 1.0f );

					UTEST( 1.0f == fValue );				
				} 
				{	
					REAL fValue = -10.f;
					Math::Clamp( fValue, 0.0f, 1.0f );

					UTEST( 0.0f == fValue );				
				} 
				{	
					REAL fValue = 0.0f;
					Math::Clamp( fValue, 0.0f, 1.0f );

					UTEST( 0.0f == fValue );				
				} 
				{	
					REAL fValue = 1.0f;
					Math::Clamp( fValue, 0.0f, 1.0f );

					UTEST( 1.0f == fValue );				
				} 
				{	
					REAL fValue = 0.5f;
					Math::Clamp( fValue, 0.0f, 1.0f );

					UTEST( 0.5f == fValue );				
				} 
			}
			// --- 3. --------------------------------------------------------------------
			LOG("-- Math::DegreeToRadian( Degree _fAngle ) && Math::RadianToDegree( Radian _fAngle ) \n");
			{
				{	
					UTEST( 0.0f				== Math::DegreeToRadian( 0.0f ) ); 
					UTEST( Math::PI / 4.0f	== Math::DegreeToRadian( 45.0f ) );	
					UTEST( Math::PI / 2.0f	== Math::DegreeToRadian( 90.0f ) );	
					UTEST( Math::PI			== Math::DegreeToRadian( 180.0f ) ); 
					UTEST( 2 * Math::PI 	== Math::DegreeToRadian( 360.0f ) );					
				}  
				{	
					UTEST( 0.0f		== Math::RadianToDegree( 0.0f ) );	
					UTEST( 45.0f	== Math::RadianToDegree( Math::PI / 4.0f ) );	
					UTEST( 90.0f	== Math::RadianToDegree( Math::PI / 2.0f ) );	 
					UTEST( 180.0f	== Math::RadianToDegree( Math::PI ) );	 
					UTEST( 360.0f	== Math::RadianToDegree( 2 * Math::PI ) );			
				}  
				{  
					UTEST( Math::Equal( -23.0f	, Math::RadianToDegree( Math::DegreeToRadian( -23.0f ) ) ) );
					UTEST( Math::Equal( 500.0f 	, Math::RadianToDegree( Math::DegreeToRadian( 500.0f ) ) ) );
					UTEST( Math::Equal( -180.0f	, Math::RadianToDegree( Math::DegreeToRadian( -180.0f ) ) ) );
				}
			}
			// --- 4. --------------------------------------------------------------------
			LOG("-- Math::Sin( Radian _fAngle ) && Math::ASin( REAL _fValue ) \n");
			{
				{	
					UTEST( Math::Equal( 0.0f, Math::Sin( 0.0f ) ) );  	
					UTEST( Math::Equal( 0.0f, Math::Sin( Math::PI ) ) );  	
					UTEST( Math::Equal( 1.0f, Math::Sin( Math::PI / 2 ) ) );  	
					UTEST( Math::Equal( 1.0f, Math::Sin( - 3 * Math::PI / 2 ) ) );  	

					UTEST( Math::Equal( 0.0f, 	Math::Sin( Math::PI * 2 ) ) );  	
					UTEST( Math::Equal( 0.0f, 	Math::Sin( Math::PI * 3 ) ) );  	
					UTEST( Math::Equal( -1.0f, 	Math::Sin( - Math::PI / 2 ) ) );  	
					UTEST( Math::Equal( -1.0f, 	Math::Sin( 3 * Math::PI / 2 ) ) );  	
				}   
				{	
					UTEST( Math::Equal( 0.0f				, Math::ASin( 0.0f ) ) );
					UTEST( Math::Equal( ( Math::PI / 2 )	, Math::ASin( 1.0f  ) ) );
					UTEST( Math::Equal( ( - Math::PI / 2 )	, Math::ASin( -1.0f ) ) );
				}   
				{	
					UTEST( Math::Equal( 0.0f				, Math::ASin( Math::Sin( 0.0f ) ) ) );  
					UTEST( Math::Equal( ( Math::PI / 4 )	, Math::ASin( Math::Sin( Math::PI / 4 ) ) ) );  
					UTEST( Math::Equal( ( Math::PI / 2 )	, Math::ASin( Math::Sin( Math::PI / 2 ) ) ) );   
					UTEST( Math::Equal( ( - Math::PI / 4 )	, Math::ASin( Math::Sin( - Math::PI / 4 ) ) ) );    
					UTEST( Math::Equal( ( - Math::PI / 2 )	, Math::ASin( Math::Sin( - Math::PI / 2 ) ) ) );   
				}   
			}
			// --- 5. --------------------------------------------------------------------
			LOG("-- Math::Cos( Radian _fAngle ) && Math::ACos( REAL _fValue ) \n");
			{
				{	
					UTEST( Math::Equal( 1.0f	, Math::Cos( 0.0f ) ) );  	
					UTEST( Math::Equal( -1.0f	, Math::Cos( Math::PI ) ) );  	
					UTEST( Math::Equal( 0.0f	, Math::Cos( Math::PI / 2 ) ) );  	
					UTEST( Math::Equal( 0.0f	, Math::Cos( - 3 * Math::PI / 2 ) ) );  	

					UTEST( Math::Equal( 1.0f	, Math::Cos( Math::PI * 2 ) ) );  	
					UTEST( Math::Equal( -1.0f	, Math::Cos( Math::PI * 3 ) ) );  	
					UTEST( Math::Equal( 0.0f	, Math::Cos( - Math::PI / 2 ) ) );  	
					UTEST( Math::Equal( 0.0f	, Math::Cos( 3 * Math::PI / 2 ) ) );
				}
				{	
					UTEST( Math::Equal( ( Math::PI / 2 ), Math::ACos( 0.0f ) ) );
					UTEST( Math::Equal( 0.0f			, Math::ACos( 1.0f  ) ) );
					UTEST( Math::Equal( Math::PI		, Math::ACos( -1.0f ) ) );
				}   
				{	
					UTEST( Math::Equal( 0.0f					, Math::ACos( Math::Cos( 0.0f ) ) ) );  
					UTEST( Math::Equal( ( Math::PI / 4 )		, Math::ACos( Math::Cos( Math::PI / 4 ) ) ) );   
					UTEST( Math::Equal( ( Math::PI / 2 )		, Math::ACos( Math::Cos( Math::PI / 2 ) ) ) );   
					UTEST( Math::Equal( ( 3 * Math::PI / 4 )	, Math::ACos( Math::Cos( 3 * Math::PI / 4 ) ) ) );   
					UTEST( Math::Equal( Math::PI				, Math::ACos( Math::Cos( Math::PI ) ) ) );   
				}   
			}
			 
			// --- 6. --------------------------------------------------------------------
			LOG("-- Math::Tan( Radian _fAngle ) && Math::ATan( REAL _fValue ) && Math::ATan2( REAL _fValueY, REAL _fValueX ) \n");
			{
				{	
					UTEST( Math::Equal( 0.0f						, Math::Tan( 0.0f ) ) );  	
					UTEST( Math::Equal( ( Math::SquareRoot( 3 )/3 )	, Math::Tan( Math::PI / 6 ) ) );  	
					UTEST( Math::Equal( 1.0f						, Math::Tan( Math::PI / 4 ) ) );  	
					UTEST( Math::Equal( Math::SquareRoot( 3 )		, Math::Tan( Math::PI / 3 ) ) );  
					
					UTEST( Math::Equal( -( Math::SquareRoot( 3 )/3 ), Math::Tan( - Math::PI / 6 ) ) );  	
					UTEST( Math::Equal( -1.0f						, Math::Tan( - Math::PI / 4 ) ) );  	
					UTEST( Math::Equal( -Math::SquareRoot( 3 )		, Math::Tan( - Math::PI / 3 ) ) );  

					UTEST( Math::Equal( -( Math::SquareRoot( 3 )/3 ), Math::Tan( Math::PI - Math::PI / 6 ) ) );  	
					UTEST( Math::Equal( -1.0f						, Math::Tan( Math::PI - Math::PI / 4 ) ) );  	
					UTEST( Math::Equal( -Math::SquareRoot( 3 )		, Math::Tan( Math::PI - Math::PI / 3 ) ) );  

					UTEST( Math::Equal( ( Math::SquareRoot( 3 )/3 )	, Math::Tan( - Math::PI + Math::PI / 6 ) ) );  	
					UTEST( Math::Equal( 1.0f						, Math::Tan( - Math::PI + Math::PI / 4 ) ) );  	
					UTEST( Math::Equal( Math::SquareRoot( 3 )		, Math::Tan( - Math::PI + Math::PI / 3 ) ) );  
				}
				{	
					UTEST( Math::Equal( 0.0f				, Math::ATan( 0.0f ) ) );
					UTEST( Math::Equal( ( Math::PI / 6 )	, Math::ATan( Math::SquareRoot( 3 )/3 ) ) );
					UTEST( Math::Equal( ( Math::PI / 4 )	, Math::ATan( 1.0f ) ) );
					UTEST( Math::Equal( ( Math::PI / 3 )	, Math::ATan( Math::SquareRoot( 3 ) ) ) );
					
					UTEST( Math::Equal( 1.373401f	, Math::ATan( 5.0f ) ) );
				}   
				{	
					UTEST( Math::Equal( 0.099669f	, Math::ATan2( 0.5f, 5.0f ) ) );
				}    
			}
			// --- 7. --------------------------------------------------------------------
			LOG("-- Math::Square( REAL ) \n");
			{
				UTEST( Math::Equal( 0.0f	, Math::Square( 0.0f ) ) );
				UTEST( Math::Equal( 0.25f	, Math::Square( 0.5f ) ) );
				UTEST( Math::Equal( 0.25f	, Math::Square( -0.5f ) ) );
				UTEST( Math::Equal( 1.0f	, Math::Square( 1.0f ) ) );
				UTEST( Math::Equal( 1.0f	, Math::Square( -1.0f ) ) );
				UTEST( Math::Equal( 4.0f	, Math::Square( 2.0f ) ) );
				UTEST( Math::Equal( 4.0f	, Math::Square( -2.0f ) ) );
			}
			// --- 8. --------------------------------------------------------------------
			LOG("-- Math::SquareRoot( REAL ) \n");
			{
				UTEST( Math::Equal( 0.0f	, Math::SquareRoot( 0.0f ) ) );
				UTEST( Math::Equal( 0.5f	, Math::SquareRoot( 0.25f ) ) );
				UTEST( Math::Equal( 1.0f	, Math::SquareRoot( 1.0f ) ) );
				UTEST( Math::Equal( 2.0f	, Math::SquareRoot( 4.0f ) ) );
			}
			// --- 9. --------------------------------------------------------------------
			LOG("-- Math::InverseSquareRoot( REAL ) \n");
			{ 
				UTEST( Math::Equal( 2.0f	, Math::InverseSquareRoot( 0.25f ) ) );
				UTEST( Math::Equal( 1.0f	, Math::InverseSquareRoot( 1.0f ) ) );
				UTEST( Math::Equal( 0.5f	, Math::InverseSquareRoot( 4.0f ) ) );
			}
			// --- 10. --------------------------------------------------------------------
			LOG("-- Math::NextMultiple( UINT32 _ulValue, UINT32 _ulMultiple ) \n");
			{
				UTEST( Math::Equal( 2.0f	, Math::NextMultiple( 1.0f, 2.0f ) ) );
				UTEST( Math::Equal( 6.0f	, Math::NextMultiple( 5.0f, 2.0f ) ) );
				UTEST( Math::Equal( 32.0f	, Math::NextMultiple( 32.0f, 32.0f ) ) );
				UTEST( Math::Equal( 64.0f	, Math::NextMultiple( 48.0f, 32.0f ) ) );
			}
			 
	 		return UTEST_RESULT; 
		}
	}
}

namespace Melting
{ 
	namespace Math
	{
		
		template< class Component, UINT32 Dim >
		bool Vector< Component, Dim >::UnitTest( void )
		{
			LOG("::::::::: CLASS : Vector< Component, Dim > ::::::::::\n");
			UTEST_START;
			
			return UTEST_RESULT;
		}

		template<>
		bool Vector2< INT16 >::UnitTest( void )
		{
			LOG("::::::::: CLASS : Vector2< INT16 > ::::::::::\n");
			UTEST_START;
			
			typedef Vector2< INT16 > V2;
			
			// operator+
			UTEST( V2(  5,  6 ) + V2(  7,  9 ) == V2( 12, 15 ) );
			UTEST( V2(  5,  6 ) + V2( -7,  9 ) == V2( -2, 15 ) );
			UTEST( V2(  5,  6 ) + V2(  7, -9 ) == V2( 12, -3 ) );
			UTEST( V2(  5,  6 ) + V2(  0,  0 ) == V2(  5,  6 ) );
			UTEST( V2( -5, -6 ) + V2(  0,  0 ) == V2( -5, -6 ) );
			
			UTEST( V2(  7,  9 ) + V2(  5,  6 ) == V2( 12, 15 ) );
			UTEST( V2( -7,  9 ) + V2(  5,  6 ) == V2( -2, 15 ) );
			UTEST( V2(  7, -9 ) + V2(  5,  6 ) == V2( 12, -3 ) );
			UTEST( V2(  0,  0 ) + V2(  5,  6 ) == V2(  5,  6 ) );
			UTEST( V2(  0,  0 ) + V2( -5, -6 ) == V2( -5, -6 ) );
			
			// operator-
			UTEST( V2(  5,  6 ) - V2(  7,  9 ) == V2( -2, -3 ) );
			UTEST( V2(  5,  6 ) - V2( -7,  9 ) == V2( 12, -3 ) );
			UTEST( V2(  5,  6 ) - V2(  7, -9 ) == V2( -2, 15 ) );
			UTEST( V2(  5,  6 ) - V2(  0,  0 ) == V2(  5,  6 ) );
			UTEST( V2( -5, -6 ) - V2(  0,  0 ) == V2( -5, -6 ) );
			
			// operator*
			UTEST( V2(  5,  6 ) * V2(  7,  9 ) == V2( 35, 54 ) );
			UTEST( V2(  5,  6 ) * V2( -7,  9 ) == V2( -35, 54 ) );
			UTEST( V2(  5,  6 ) * V2(  7, -9 ) == V2( 35, -54 ) );
			UTEST( V2(  0,  0 ) * V2(  5,  6 ) == V2(  0,  0 ) );
			UTEST( V2( -5, -6 ) * V2(  0,  0 ) == V2(  0,  0 ) );
			
			// operator/
			UTEST( V2(  5,  6 ) / V2(  7,  9 ) == V2(  0,  0 ) );
			UTEST( V2(  7,  9 ) / V2(  5,  6 ) == V2(  1,  1 ) );
			UTEST( V2( 15,  9 ) / V2(  5,  6 ) == V2(  3,  1 ) );
			UTEST( V2(  5,  6 ) / V2( -7,  9 ) == V2(  0,  0 ) );
			UTEST( V2( 43, 24 ) / V2( -7,  9 ) == V2( -6,  2 ) );
			
			return UTEST_RESULT;
		}
		
		template<>
		bool Vector2< REAL >::UnitTest( void )
		{
			LOG("::::::::: CLASS : Vector2< REAL > ::::::::::\n");
			UTEST_START;
			
			typedef Vector2< REAL > V2;
			typedef REAL Fx;
			
			// operator+
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) + V2( Fx( 7 ), Fx( 9 ) ) == V2( Fx(12 ), Fx(15 ) ) );
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) + V2( Fx(-7 ), Fx( 9 ) ) == V2( Fx(-2 ), Fx(15 ) ) );
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) + V2( Fx( 7 ), Fx(-9 ) ) == V2( Fx(12 ), Fx(-3 ) ) );
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) + V2( Fx( 0 ), Fx( 0 ) ) == V2( Fx( 5 ), Fx( 6 ) ) );
			UTEST( V2( Fx(-5 ), Fx(-6 ) ) + V2( Fx( 0 ), Fx( 0 ) ) == V2( Fx(-5 ), Fx(-6 ) ) );
			
			UTEST( V2( Fx( 7 ), Fx( 9 ) ) + V2( Fx( 5 ), Fx( 6 ) ) == V2( Fx(12 ), Fx(15 ) ) );
			UTEST( V2( Fx(-7 ), Fx( 9 ) ) + V2( Fx( 5 ), Fx( 6 ) ) == V2( Fx(-2 ), Fx(15 ) ) );
			UTEST( V2( Fx( 7 ), Fx(-9 ) ) + V2( Fx( 5 ), Fx( 6 ) ) == V2( Fx(12 ), Fx(-3 ) ) );
			UTEST( V2( Fx( 0 ), Fx( 0 ) ) + V2( Fx( 5 ), Fx( 6 ) ) == V2( Fx( 5 ), Fx( 6 ) ) );
			UTEST( V2( Fx( 0 ), Fx( 0 ) ) + V2( Fx(-5 ), Fx(-6 ) ) == V2( Fx(-5 ), Fx(-6 ) ) );
			
			// operator-
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) - V2( Fx( 7 ), Fx( 9 ) ) == V2( Fx(-2 ), Fx(-3 ) ) );
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) - V2( Fx(-7 ), Fx( 9 ) ) == V2( Fx(12 ), Fx(-3 ) ) );
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) - V2( Fx( 7 ), Fx(-9 ) ) == V2( Fx(-2 ), Fx(15 ) ) );
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) - V2( Fx( 0 ), Fx( 0 ) ) == V2( Fx( 5 ), Fx( 6 ) ) );
			UTEST( V2( Fx(-5 ), Fx(-6 ) ) - V2( Fx( 0 ), Fx( 0 ) ) == V2( Fx(-5 ), Fx(-6 ) ) );
			
			// operator*
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) * V2( Fx( 7 ), Fx( 9 ) ) == V2( Fx(35 ), Fx(54 ) ) );
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) * V2( Fx(-7 ), Fx( 9 ) ) == V2( Fx(-35), Fx(54 ) ) );
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) * V2( Fx( 7 ), Fx(-9 ) ) == V2( Fx(35 ), Fx(-54) ) );
			UTEST( V2( Fx( 0 ), Fx( 0 ) ) * V2( Fx( 5 ), Fx( 6 ) ) == V2( Fx( 0 ), Fx( 0 ) ) );
			UTEST( V2( Fx(-5 ), Fx(-6 ) ) * V2( Fx( 0 ), Fx( 0 ) ) == V2( Fx( 0 ), Fx( 0 ) ) );
			
			// operator/
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) / V2( Fx( 7 ), Fx( 9 ) ) == V2( Fx(  5 ) / Fx(  7 ),  Fx(  6 ) / Fx( 9 ) ) );
			UTEST( V2( Fx( 7 ), Fx( 9 ) ) / V2( Fx( 5 ), Fx( 6 ) ) == V2( Fx(  7 ) / Fx(  5 ),  Fx(  9 ) / Fx( 6 ) ) );
			UTEST( V2( Fx(15 ), Fx( 9 ) ) / V2( Fx( 5 ), Fx( 6 ) ) == V2( Fx( 15 ) / Fx(  5 ),  Fx(  9 ) / Fx( 6 ) ) );
			UTEST( V2( Fx( 5 ), Fx( 6 ) ) / V2( Fx(-7 ), Fx( 9 ) ) == V2( Fx(  5 ) / Fx( -7 ),  Fx(  6 ) / Fx( 9 ) ) );
			UTEST( V2( Fx(43 ), Fx(24 ) ) / V2( Fx(-7 ), Fx( 9 ) ) == V2( Fx( 43 ) / Fx( -7 ),  Fx( 24 ) / Fx( 9 ) ) );
			
			return UTEST_RESULT;
		}
	}// namespace Math
}// namespace Melting



namespace Melting
{
	namespace Math
	{
		bool Mtx_UnitTest( void )
		{
			LOG( "::::::::: Functions : Matrix  ::::::::::\n" );

			UTEST_START; 
			
			// --- 1. --------------------------------------------------------------------
			LOG( "-- Mtx_BuildNull( )  \n" );
			{	
				Matrix44 mtx;

				Mtx_BuildNull( mtx );

				Matrix44 mtxExpected;
				mtxExpected.m[0][0] = 0.0f;
				mtxExpected.m[0][1] = 0.0f;
				mtxExpected.m[0][2] = 0.0f;
				mtxExpected.m[0][3] = 0.0f;

				mtxExpected.m[1][0] = 0.0f;
				mtxExpected.m[1][1] = 0.0f;
				mtxExpected.m[1][2] = 0.0f;
				mtxExpected.m[1][3] = 0.0f;

				mtxExpected.m[2][0] = 0.0f;
				mtxExpected.m[2][1] = 0.0f;
				mtxExpected.m[2][2] = 0.0f;
				mtxExpected.m[2][3] = 0.0f;

				mtxExpected.m[3][0] = 0.0f;
				mtxExpected.m[3][1] = 0.0f;
				mtxExpected.m[3][2] = 0.0f;
				mtxExpected.m[3][3] = 0.0f;

				UTEST( Mtx_Equal( mtx, mtxExpected ) );
			}
			
			// --- 1. --------------------------------------------------------------------
			LOG( "-- Mtx_BuildIdentity( )  \n" );
			{	
				Matrix44 mtx;

				Mtx_BuildIdentity( mtx );

				Matrix44 mtxExpected;
				mtxExpected.m[0][0] = 1.0f;
				mtxExpected.m[0][1] = 0.0f;
				mtxExpected.m[0][2] = 0.0f;
				mtxExpected.m[0][3] = 0.0f;

				mtxExpected.m[1][0] = 0.0f;
				mtxExpected.m[1][1] = 1.0f;
				mtxExpected.m[1][2] = 0.0f;
				mtxExpected.m[1][3] = 0.0f;

				mtxExpected.m[2][0] = 0.0f;
				mtxExpected.m[2][1] = 0.0f;
				mtxExpected.m[2][2] = 1.0f;
				mtxExpected.m[2][3] = 0.0f;

				mtxExpected.m[3][0] = 0.0f;
				mtxExpected.m[3][1] = 0.0f;
				mtxExpected.m[3][2] = 0.0f;
				mtxExpected.m[3][3] = 1.0f;

				UTEST( Mtx_Equal( mtx, mtxExpected ) );
			}
			
			// --- 2. --------------------------------------------------------------------
			LOG( "-- Mtx_BuildTranslation( )  \n" );
			{	
				Matrix44 mtx;

				Vec3 vTrans( 1.0f, 2.0f, 3.0f );

				Mtx_BuildTranslation( mtx, vTrans );

				Matrix44 mtxExpected;
				mtxExpected.m[0][0] = 1.0f;
				mtxExpected.m[0][1] = 0.0f;
				mtxExpected.m[0][2] = 0.0f;
				mtxExpected.m[0][3] = 0.0f;

				mtxExpected.m[1][0] = 0.0f;
				mtxExpected.m[1][1] = 1.0f;
				mtxExpected.m[1][2] = 0.0f;
				mtxExpected.m[1][3] = 0.0f;

				mtxExpected.m[2][0] = 0.0f;
				mtxExpected.m[2][1] = 0.0f;
				mtxExpected.m[2][2] = 1.0f;
				mtxExpected.m[2][3] = 0.0f;

				mtxExpected.m[3][0] = vTrans[0];
				mtxExpected.m[3][1] = vTrans[1];
				mtxExpected.m[3][2] = vTrans[2];
				mtxExpected.m[3][3] = 1.0f;

				UTEST( Mtx_Equal( mtx, mtxExpected ) );
			}

			// --- 3. --------------------------------------------------------------------
			LOG( "-- Mtx_BuildScale( )  \n" );
			{	
				Matrix44 mtx;

				Vec3 vScale( 1.0f, 2.0f, 3.0f );

				Mtx_BuildScale( mtx, vScale );

				Matrix44 mtxExpected;
				mtxExpected.m[0][0] = vScale[0];
				mtxExpected.m[0][1] = 0.0f;
				mtxExpected.m[0][2] = 0.0f;
				mtxExpected.m[0][3] = 0.0f;

				mtxExpected.m[1][0] = 0.0f;
				mtxExpected.m[1][1] = vScale[1];
				mtxExpected.m[1][2] = 0.0f;
				mtxExpected.m[1][3] = 0.0f;

				mtxExpected.m[2][0] = 0.0f;
				mtxExpected.m[2][1] = 0.0f;
				mtxExpected.m[2][2] = vScale[2];
				mtxExpected.m[2][3] = 0.0f;

				mtxExpected.m[3][0] = 0.0f;
				mtxExpected.m[3][1] = 0.0f;
				mtxExpected.m[3][2] = 0.0f;
				mtxExpected.m[3][3] = 1.0f;

				UTEST( Mtx_Equal( mtx, mtxExpected ) );
			}  

			// --- 4. --------------------------------------------------------------------
			LOG( "-- Mtx_Scale( )  \n" );
			{	
				Matrix44 mtxBase;
				Mtx_BuildIdentity( mtxBase ); 

				Matrix44 mtx;

				REAL fScale = 3.0f;

				Mtx_Scale( mtx, mtxBase, fScale ); 

				Matrix44 mtxExpected;
				mtxExpected.m[0][0] = fScale;
				mtxExpected.m[0][1] = 0.0f;
				mtxExpected.m[0][2] = 0.0f;
				mtxExpected.m[0][3] = 0.0f;

				mtxExpected.m[1][0] = 0.0f;
				mtxExpected.m[1][1] = fScale;
				mtxExpected.m[1][2] = 0.0f;
				mtxExpected.m[1][3] = 0.0f;

				mtxExpected.m[2][0] = 0.0f;
				mtxExpected.m[2][1] = 0.0f;
				mtxExpected.m[2][2] = fScale;
				mtxExpected.m[2][3] = 0.0f;

				mtxExpected.m[3][0] = 0.0f;
				mtxExpected.m[3][1] = 0.0f;
				mtxExpected.m[3][2] = 0.0f;
				mtxExpected.m[3][3] = fScale;

				UTEST( Mtx_Equal( mtx, mtxExpected ) );
			}   

			// --- 5. --------------------------------------------------------------------
			LOG( "-- Mtx_BuildRotationX( )  \n" );
			{
				{
					Matrix44 mtx;

					REAL fAngle = Math::PI / 2.0f;

					Mtx_BuildRotationX( mtx, fAngle );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 1.0f;
					mtxExpected.m[0][1] = 0.0f;
					mtxExpected.m[0][2] = 0.0f;
					mtxExpected.m[0][3] = 0.0f;

					mtxExpected.m[1][0] = 0.0f;
					mtxExpected.m[1][1] = 0.0f;
					mtxExpected.m[1][2] = 1.0f;
					mtxExpected.m[1][3] = 0.0f;

					mtxExpected.m[2][0] = 0.0f;
					mtxExpected.m[2][1] = -1.0f;
					mtxExpected.m[2][2] = 0.0f;
					mtxExpected.m[2][3] = 0.0f;

					mtxExpected.m[3][0] = 0.0f;
					mtxExpected.m[3][1] = 0.0f;
					mtxExpected.m[3][2] = 0.0f;
					mtxExpected.m[3][3] = 1.0f;

					UTEST( Mtx_Equal( mtx, mtxExpected ) );
				}
			}
			// --- 6. --------------------------------------------------------------------
			LOG( "-- Mtx_BuildRotationY( )  \n" );
			{
				{
					Matrix44 mtx;

					REAL fAngle = Math::PI / 2.0f;

					Mtx_BuildRotationY( mtx, fAngle );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 0.0f;
					mtxExpected.m[0][1] = 0.0f;
					mtxExpected.m[0][2] = -1.0f;
					mtxExpected.m[0][3] = 0.0f;

					mtxExpected.m[1][0] = 0.0f;
					mtxExpected.m[1][1] = 1.0f;
					mtxExpected.m[1][2] = 0.0f;
					mtxExpected.m[1][3] = 0.0f;

					mtxExpected.m[2][0] = 1.0f;
					mtxExpected.m[2][1] = 0.0f;
					mtxExpected.m[2][2] = 0.0f;
					mtxExpected.m[2][3] = 0.0f;

					mtxExpected.m[3][0] = 0.0f;
					mtxExpected.m[3][1] = 0.0f;
					mtxExpected.m[3][2] = 0.0f;
					mtxExpected.m[3][3] = 1.0f;

					UTEST( Mtx_Equal( mtx, mtxExpected ) );
				}
			}
			// --- 7. --------------------------------------------------------------------
			LOG( "-- Mtx_BuildRotationZ( )  \n" );
			{
				{
					Matrix44 mtx;

					REAL fAngle = Math::PI / 2.0f;

					Mtx_BuildRotationZ( mtx, fAngle );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 0.0f;
					mtxExpected.m[0][1] = 1.0f;
					mtxExpected.m[0][2] = 0.0f;
					mtxExpected.m[0][3] = 0.0f;

					mtxExpected.m[1][0] = -1.0f;
					mtxExpected.m[1][1] = 0.0f;
					mtxExpected.m[1][2] = 0.0f;
					mtxExpected.m[1][3] = 0.0f;

					mtxExpected.m[2][0] = 0.0f;
					mtxExpected.m[2][1] = 0.0f;
					mtxExpected.m[2][2] = 1.0f;
					mtxExpected.m[2][3] = 0.0f;

					mtxExpected.m[3][0] = 0.0f;
					mtxExpected.m[3][1] = 0.0f;
					mtxExpected.m[3][2] = 0.0f;
					mtxExpected.m[3][3] = 1.0f;

					UTEST( Mtx_Equal( mtx, mtxExpected ) );
				}
			}
			// --- 8. --------------------------------------------------------------------
			LOG( "-- Mtx_BuildRotation( )  \n" );
			{	
				{
					Matrix44 mtx;

					Vec3 vAxis	( 0.0f, 1.0f, 0.0f );
					REAL fAngle = 0;

					Mtx_BuildRotation( mtx, vAxis, fAngle );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 1.0f;
					mtxExpected.m[0][1] = 0.0f;
					mtxExpected.m[0][2] = 0.0f;
					mtxExpected.m[0][3] = 0.0f;

					mtxExpected.m[1][0] = 0.0f;
					mtxExpected.m[1][1] = 1.0f;
					mtxExpected.m[1][2] = 0.0f;
					mtxExpected.m[1][3] = 0.0f;

					mtxExpected.m[2][0] = 0.0f;
					mtxExpected.m[2][1] = 0.0f;
					mtxExpected.m[2][2] = 1.0f;
					mtxExpected.m[2][3] = 0.0f;

					mtxExpected.m[3][0] = 0.0f;
					mtxExpected.m[3][1] = 0.0f;
					mtxExpected.m[3][2] = 0.0f;
					mtxExpected.m[3][3] = 1.0f;

					UTEST( Mtx_Equal( mtx, mtxExpected ) );
				}
				{
					Matrix44 mtx;

					Vec3 vAxis	( 0.0f, 1.0f, 0.0f );
					REAL fAngle = Math::PI / 2.0f;

					Mtx_BuildRotation( mtx, vAxis, fAngle );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 0.0f;
					mtxExpected.m[0][1] = 0.0f;
					mtxExpected.m[0][2] = -1.0f;
					mtxExpected.m[0][3] = 0.0f;

					mtxExpected.m[1][0] = 0.0f;
					mtxExpected.m[1][1] = 1.0f;
					mtxExpected.m[1][2] = 0.0f;
					mtxExpected.m[1][3] = 0.0f;

					mtxExpected.m[2][0] = 1.0f;
					mtxExpected.m[2][1] = 0.0f;
					mtxExpected.m[2][2] = 0.0f;
					mtxExpected.m[2][3] = 0.0f;

					mtxExpected.m[3][0] = 0.0f;
					mtxExpected.m[3][1] = 0.0f;
					mtxExpected.m[3][2] = 0.0f;
					mtxExpected.m[3][3] = 1.0f;

					UTEST( Mtx_Equal( mtx, mtxExpected ) );
				}
				{
					Matrix44 mtx;

					Vec3 vAxis	( 0.0f, 1.0f, 0.0f );
					REAL fAngle = Math::PI;

					Mtx_BuildRotation( mtx, vAxis, fAngle );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = -1.0f;
					mtxExpected.m[0][1] = 0.0f;
					mtxExpected.m[0][2] = 0.0f;
					mtxExpected.m[0][3] = 0.0f;

					mtxExpected.m[1][0] = 0.0f;
					mtxExpected.m[1][1] = 1.0f;
					mtxExpected.m[1][2] = 0.0f;
					mtxExpected.m[1][3] = 0.0f;

					mtxExpected.m[2][0] = 0.0f;
					mtxExpected.m[2][1] = 0.0f;
					mtxExpected.m[2][2] = -1.0f;
					mtxExpected.m[2][3] = 0.0f;

					mtxExpected.m[3][0] = 0.0f;
					mtxExpected.m[3][1] = 0.0f;
					mtxExpected.m[3][2] = 0.0f;
					mtxExpected.m[3][3] = 1.0f;

					UTEST( Mtx_Equal( mtx, mtxExpected ) );
				}
				{
					Matrix44 mtx;

					Vec3 vAxis	( 0.0f, 0.0f, 1.0f );
					REAL fAngle = Math::PI;

					Mtx_BuildRotation( mtx, vAxis, fAngle );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = -1.0f;
					mtxExpected.m[0][1] = 0.0f;
					mtxExpected.m[0][2] = 0.0f;
					mtxExpected.m[0][3] = 0.0f;

					mtxExpected.m[1][0] = 0.0f;
					mtxExpected.m[1][1] = -1.0f;
					mtxExpected.m[1][2] = 0.0f;
					mtxExpected.m[1][3] = 0.0f;

					mtxExpected.m[2][0] = 0.0f;
					mtxExpected.m[2][1] = 0.0f;
					mtxExpected.m[2][2] = 1.0f;
					mtxExpected.m[2][3] = 0.0f;

					mtxExpected.m[3][0] = 0.0f;
					mtxExpected.m[3][1] = 0.0f;
					mtxExpected.m[3][2] = 0.0f;
					mtxExpected.m[3][3] = 1.0f;

					UTEST( Mtx_Equal( mtx, mtxExpected ) );
				}
				{
					Matrix44 mtx;

					Vec3 vAxis	(  1.0f, 0.0f, 0.0f );
					REAL fAngle = Math::PI / 2.0f;

					Mtx_BuildRotation( mtx, vAxis, fAngle );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 1.0f;
					mtxExpected.m[0][1] = 0.0f;
					mtxExpected.m[0][2] = 0.0f;
					mtxExpected.m[0][3] = 0.0f;

					mtxExpected.m[1][0] = 0.0f;
					mtxExpected.m[1][1] = 0.0f;
					mtxExpected.m[1][2] = 1.0f;
					mtxExpected.m[1][3] = 0.0f;

					mtxExpected.m[2][0] = 0.0f;
					mtxExpected.m[2][1] = -1.0f;
					mtxExpected.m[2][2] = 0.0f;
					mtxExpected.m[2][3] = 0.0f;

					mtxExpected.m[3][0] = 0.0f;
					mtxExpected.m[3][1] = 0.0f;
					mtxExpected.m[3][2] = 0.0f;
					mtxExpected.m[3][3] = 1.0f;

					UTEST( Mtx_Equal( mtx, mtxExpected ) );
				}
				{
					Matrix44 mtx;

					Vec3 vAxis	(  1.0f, 0.0f, 0.0f );
					REAL fAngle = Math::PI;

					Mtx_BuildRotation( mtx, vAxis, fAngle );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 1.0f;
					mtxExpected.m[0][1] = 0.0f;
					mtxExpected.m[0][2] = 0.0f;
					mtxExpected.m[0][3] = 0.0f;

					mtxExpected.m[1][0] = 0.0f;
					mtxExpected.m[1][1] = -1.0f;
					mtxExpected.m[1][2] = 0.0f;
					mtxExpected.m[1][3] = 0.0f;

					mtxExpected.m[2][0] = 0.0f;
					mtxExpected.m[2][1] = 0.0f;
					mtxExpected.m[2][2] = -1.0f;
					mtxExpected.m[2][3] = 0.0f;

					mtxExpected.m[3][0] = 0.0f;
					mtxExpected.m[3][1] = 0.0f;
					mtxExpected.m[3][2] = 0.0f;
					mtxExpected.m[3][3] = 1.0f;

					UTEST( Mtx_Equal( mtx, mtxExpected ) );
				}
				 
				{
					Matrix44 mtxA;
					Mtx_BuildIdentity( mtxA );
					Matrix44 mtxB;
					Mtx_BuildIdentity( mtxB );

					// 90° around X
					{
						Matrix44 mtxRotation;
						Vec3 vAxis	( 1.0f, 0.0f, 0.0f );
						REAL fAngle = 90.f * Math::PI/180.0f;

						//Mtx_BuildRotation( mtxRotation, vAxis, fAngle );
						Mtx_BuildRotationX( mtxRotation, fAngle );

						Mtx_Concat( mtxA, mtxRotation, mtxB );
					}
					// 40° around Y
					{
						Matrix44 mtxRotation;
						Vec3 vAxis	( 0.0f, 1.0f, 0.0f );
						REAL fAngle = 40.f * Math::PI/180.0f;

						//Mtx_BuildRotation( mtxRotation, vAxis, fAngle );
						Mtx_BuildRotationY( mtxRotation, fAngle );

						Mtx_Concat( mtxB , mtxRotation, mtxA );
					}
					// 50° around Z
					{
						Matrix44 mtxRotation;
						Vec3 vAxis	( 0.0f, 0.0f, 1.0f );
						REAL fAngle = 50.f * Math::PI/180.0f;

						//Mtx_BuildRotation( mtxRotation, vAxis, fAngle );
						Mtx_BuildRotationZ( mtxRotation, fAngle );

						Mtx_Concat( mtxA, mtxRotation, mtxB );
					}

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 0.492403876506104f;
					mtxExpected.m[0][1] = 0.586824088833465f;
					mtxExpected.m[0][2] = -0.642787609686539f;
					mtxExpected.m[0][3] = 0.0f;

					mtxExpected.m[1][0] = 0.413175911166535f;
					mtxExpected.m[1][1] = 0.492403876506104f;
					mtxExpected.m[1][2] = 0.766044443118978f;
					mtxExpected.m[1][3] = 0.0f;

					mtxExpected.m[2][0] = 0.766044443118978f;
					mtxExpected.m[2][1] = -0.642787609686539f;
					mtxExpected.m[2][2] = 0.0f;
					mtxExpected.m[2][3] = 0.0f; 

					mtxExpected.m[3][0] = 0.0f;
					mtxExpected.m[3][1] = 0.0f;
					mtxExpected.m[3][2] = 0.0f;
					mtxExpected.m[3][3] = 1.0f;

					UTEST( Mtx_Equal( mtxB, mtxExpected ) );
				}
			} 
			
			// --- 9. --------------------------------------------------------------------
			LOG("-- Mtx_Transpose( ) \n");
			{	
				{
					Matrix44 mtx;
					mtx.m[0][0] = 0;
					mtx.m[0][1] = 1;
					mtx.m[0][2] = 2;
					mtx.m[0][3] = 3;

					mtx.m[1][0] = 4;
					mtx.m[1][1] = 5;
					mtx.m[1][2] = 6;
					mtx.m[1][3] = 7;

					mtx.m[2][0] = 8;
					mtx.m[2][1] = 9;
					mtx.m[2][2] = 10;
					mtx.m[2][3] = 11; 

					mtx.m[3][0] = 12;
					mtx.m[3][1] = 13;
					mtx.m[3][2] = 14;
					mtx.m[3][3] = 15;
					
					Matrix44 mtxResult;
					Mtx_Transpose( mtx, mtxResult );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 0;
					mtxExpected.m[0][1] = 4;
					mtxExpected.m[0][2] = 8;
					mtxExpected.m[0][3] = 12;

					mtxExpected.m[1][0] = 1;
					mtxExpected.m[1][1] = 5;
					mtxExpected.m[1][2] = 9;
					mtxExpected.m[1][3] = 13;

					mtxExpected.m[2][0] = 2;
					mtxExpected.m[2][1] = 6;
					mtxExpected.m[2][2] = 10;
					mtxExpected.m[2][3] = 14; 

					mtxExpected.m[3][0] = 3;
					mtxExpected.m[3][1] = 7;
					mtxExpected.m[3][2] = 11;
					mtxExpected.m[3][3] = 15;
					
					UTEST( Mtx_Equal( mtxResult, mtxExpected ) );
				}
			}
			
			// --- 10. --------------------------------------------------------------------
			LOG("-- Mtx_Inverse( ) \n");
			{	
				// inverse of a rotation matrix is its transpose, because the matrix is orthogonal
				{
					Matrix44 mtx;
					Mtx_BuildRotationX( mtx, Math::PI / 2.0f );

					Matrix44 mtxTrans;
					Mtx_Transpose( mtx, mtxTrans );

					Vec3 vDet;
					Matrix44 mtxInv;
					Mtx_Inverse( mtx, mtxInv, vDet );

					UTEST( Mtx_Equal( mtxInv, mtxTrans ) );
				}
				// inverse of a translation matrix is the same matrix with the opposite translation
				{
					Matrix44 mtx;
					Vec3 vT ( 1, 2, 3 );
					Mtx_BuildTranslation( mtx, vT );

					Matrix44 mtx2;
					Vec3 vT2 ( -1, -2, -3 );
					Mtx_BuildTranslation( mtx2, vT2 );

					Vec3 vDet;
					Matrix44 mtxInv;
					Mtx_Inverse( mtx, mtxInv, vDet );

					UTEST( Mtx_Equal( mtxInv, mtx2 ) );
				}
			}
			
			// --- 11. --------------------------------------------------------------------
			LOG("-- Mtx_Add( ) \n");
			{	
				{
					Matrix44 mtx;
					mtx.m[0][0] = 0;
					mtx.m[0][1] = 1;
					mtx.m[0][2] = 2;
					mtx.m[0][3] = 3;

					mtx.m[1][0] = 4;
					mtx.m[1][1] = 5;
					mtx.m[1][2] = 6;
					mtx.m[1][3] = 7;

					mtx.m[2][0] = 8;
					mtx.m[2][1] = 9;
					mtx.m[2][2] = 10;
					mtx.m[2][3] = 11; 

					mtx.m[3][0] = 12;
					mtx.m[3][1] = 13;
					mtx.m[3][2] = 14;
					mtx.m[3][3] = 15;
					
					Matrix44 mtxResult;
					Mtx_Add( mtx, mtx, mtxResult );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 0;
					mtxExpected.m[0][1] = 2;
					mtxExpected.m[0][2] = 4;
					mtxExpected.m[0][3] = 6;

					mtxExpected.m[1][0] = 8;
					mtxExpected.m[1][1] = 10;
					mtxExpected.m[1][2] = 12;
					mtxExpected.m[1][3] = 14;

					mtxExpected.m[2][0] = 16;
					mtxExpected.m[2][1] = 18;
					mtxExpected.m[2][2] = 20;
					mtxExpected.m[2][3] = 22; 

					mtxExpected.m[3][0] = 24;
					mtxExpected.m[3][1] = 26;
					mtxExpected.m[3][2] = 28;
					mtxExpected.m[3][3] = 30;
					
					UTEST( Mtx_Equal( mtxResult, mtxExpected ) );
				}
			}
			
			// --- 12. --------------------------------------------------------------------
			LOG("-- Mtx_Concat( ) \n");
			{	
				{
					Matrix44 mtx;
					Mtx_BuildIdentity( mtx );
					
					Matrix44 mtxResult;
					Mtx_Concat( mtx, mtx, mtxResult );

					Matrix44 mtxExpected;
					Mtx_BuildIdentity( mtxExpected ); 
					
					UTEST( Mtx_Equal( mtxResult, mtxExpected ) );
				}
				{
					Matrix44 mtx;
					mtx.m[0][0] = 0;
					mtx.m[0][1] = 1;
					mtx.m[0][2] = 2;
					mtx.m[0][3] = 3;

					mtx.m[1][0] = 4;
					mtx.m[1][1] = 5;
					mtx.m[1][2] = 6;
					mtx.m[1][3] = 7;

					mtx.m[2][0] = 8;
					mtx.m[2][1] = 9;
					mtx.m[2][2] = 10;
					mtx.m[2][3] = 11; 

					mtx.m[3][0] = 12;
					mtx.m[3][1] = 13;
					mtx.m[3][2] = 14;
					mtx.m[3][3] = 15;
					
					Matrix44 mtxResult;
					Mtx_Concat( mtx, mtx, mtxResult );

					Matrix44 mtxExpected;
					mtxExpected.m[0][0] = 56;
					mtxExpected.m[0][1] = 62;
					mtxExpected.m[0][2] = 68;
					mtxExpected.m[0][3] = 74;

					mtxExpected.m[1][0] = 152;
					mtxExpected.m[1][1] = 174;
					mtxExpected.m[1][2] = 196;
					mtxExpected.m[1][3] = 218;

					mtxExpected.m[2][0] = 248;
					mtxExpected.m[2][1] = 286;
					mtxExpected.m[2][2] = 324;
					mtxExpected.m[2][3] = 362; 

					mtxExpected.m[3][0] = 344;
					mtxExpected.m[3][1] = 398;
					mtxExpected.m[3][2] = 452;
					mtxExpected.m[3][3] = 506;
					
					UTEST( Mtx_Equal( mtxResult, mtxExpected ) );
				}
			} 
			
			// --- 13. --------------------------------------------------------------------
			LOG("-- Mtx_BuildFromSRT( ) \n");
			{	 
				Vec3 vT( 1,2,3 );
				Quaternion qR( Math::PI / 2.0f, Vec3( 0, 1, 0 ) );
				Vec3 vS( 2,2,2 );

				Matrix44 mtxResult;
				Mtx_BuildFromSRT( mtxResult, vT, qR, vS );

				Matrix44 mtxExpected;
				mtxExpected.m[0][0] = 0;
				mtxExpected.m[0][1] = 0;
				mtxExpected.m[0][2] = -2;
				mtxExpected.m[0][3] = 0;

				mtxExpected.m[1][0] = 0;
				mtxExpected.m[1][1] = 2;
				mtxExpected.m[1][2] = 0;
				mtxExpected.m[1][3] = 0;

				mtxExpected.m[2][0] = 2;
				mtxExpected.m[2][1] = 0;
				mtxExpected.m[2][2] = 0;
				mtxExpected.m[2][3] = 0; 

				mtxExpected.m[3][0] = vT[0];
				mtxExpected.m[3][1] = vT[1];
				mtxExpected.m[3][2] = vT[2];
				mtxExpected.m[3][3] = 1;

				UTEST( Mtx_Equal( mtxResult, mtxExpected ) ); 
			} 
			// --- 14. --------------------------------------------------------------------
			LOG("-- Mtx_BuildPerspectiveRH( ) - TODO \n");
			{	 
			} 
			// --- 15. --------------------------------------------------------------------
			LOG("-- Mtx_BuildOrtographicRH( ) - TODO \n");
			{	 
			} 
			// --- 16. --------------------------------------------------------------------
			LOG("-- Mtx_MultVector( ) \n");
			{	 
				Vec3 vOriginal( 0,0,0 );

				Matrix44 mtx;
				Mtx_BuildTranslation( mtx, Vec3(10,1,-1) );

				Vec3 vResult;
				Mtx_MultVector( mtx, vOriginal, vResult );

				Vec3 vExpected; 

				UTEST( vResult == Vec3(10,1,-1) ); 
			}  
 
			return UTEST_RESULT;
		}
 
		bool SIMDUnitTest( void )
		{
			LOG("::::::::: Functions : SIMD ::::::::::\n");

			UTEST_START;
			{
				//float tmp[4] = {3,2,1,0};

				// // Method 1
				//SIMDVector v = _mm_set_ps(tmp[3],tmp[2],tmp[1],tmp[0]);

				// // Method 2
				//SIMDVector u = _mm_loadu_ps( tmp );
				//u = _mm_shuffle_ps( u, u, _MM_SHUFFLE(2,2,1,0) );
			}
			
			// --- 1. --------------------------------------------------------------------
			LOG("-- Load(  ) && Store( )  \n");
			{	
				{
					Vec3 vOriginal( 0, 1, 2 );
					SIMDVector input;

					SIMD_LOAD_A_FLOAT3( vOriginal, input ); 


					Vec3 vResult;
					
					SIMD_STORE_A_FLOAT3( vResult, input ); 
					
					UTEST( vResult == vOriginal );					
				}
				{
					float vOriginal[4] = { 0, 1, 2, 3 };
					SIMDVector input;

					SIMD_LOAD_A_FLOAT4( vOriginal, input ); 


					ALIGN_16( float vResult[4]);
					
					SIMD_STORE_A_FLOAT4( vResult, input ); 
					
					UTEST( vResult[0] == vOriginal[0] );					
					UTEST( vResult[1] == vOriginal[1] );					
					UTEST( vResult[2] == vOriginal[2] );					
					UTEST( vResult[3] == vOriginal[3] );					
				}
				{
					float vOriginal[4] = { 0, 1, 2, 3 };
					SIMDVector input;

					SIMD_LOAD_A_4FLOAT( vOriginal[0], vOriginal[1], vOriginal[2], vOriginal[3] , input ); 


					ALIGN_16( float vResult[4] );
					
					SIMD_STORE_A_FLOAT4( vResult, input ); 
					
					UTEST( vResult[0] == vOriginal[0] );					
					UTEST( vResult[1] == vOriginal[1] );					
					UTEST( vResult[2] == vOriginal[2] );					
					UTEST( vResult[3] == vOriginal[3] );										
				}
			}				
			
			// --- 2. --------------------------------------------------------------------
			LOG("-- Add(  ) \n");
			{	
				Vec3 vA( 0, 1, 2 );
				Vec3 vB( 9, 10, 20 );
				SIMDVector a, b;

				SIMD_LOAD_A_FLOAT3( vA, a );
				SIMD_LOAD_A_FLOAT3( vB, b );

				SIMDVector result = _mm_add_ps( a, b );

				Vec3 vResult;
				
				SIMD_STORE_A_FLOAT3( vResult, result ); 
				
				UTEST( vResult == Vec3( 9, 11, 22 ) );
			}
			
			// --- 3. --------------------------------------------------------------------
			LOG("-- DotProduct( ) \n");
			{	 
				Vec3 vA( 0, 1, 2 );
				Vec3 vB( 3, 4, 5 );
				SIMDVector a, b;
				SIMD_LOAD_A_FLOAT3_AS_0ZYX( vA, a );
				SIMD_LOAD_A_FLOAT3_AS_0ZYX( vB, b );
				
				SIMDVector result;
				SIMD_DOT_PRODUCT( a, b, result );

				Vec3 vResult;				
				SIMD_STORE_A_FLOAT3( vResult, result ); 
				
				UTEST( vResult[0] == 14 );
			}
			// --- 4. --------------------------------------------------------------------
			LOG("-- CrossProduct( ) \n");
			{	 		
				Vec3 vA( 0, 1, 2 );
				Vec3 vB( 3, 4, 5 );
				SIMDVector a, b;
				SIMD_LOAD_A_FLOAT3( vA, a );
				SIMD_LOAD_A_FLOAT3( vB, b );
				  
				SIMDVector result;
				SIMD_CROSS_PRODUCT_3( a, b, result );

				Vec3 vResult;				
				SIMD_STORE_A_FLOAT3( vResult, result ); 
				
				UTEST( vResult == Vec3( -3, 6, -3 ) );
			}
			
			// --- 5. --------------------------------------------------------------------
			LOG("-- SqrLength( ) \n");
			{	 		
				Vec3 vA ( 3, 4, 5 );
				SIMDVector a;
				SIMD_LOAD_A_FLOAT3_AS_0ZYX( vA, a ); 
				  
				SIMDVector result;
				SIMD_SQRLENGTH( a, result );

				Vec3 vResult;				
				SIMD_STORE_A_FLOAT3( vResult, result ); 
				
				UTEST( vResult[0] == 50 );
			}

			// --- 6. --------------------------------------------------------------------
			LOG("-- Length( ) \n");
			{	 		
				Vec3 vA ( 3, 4, 5 );
				SIMDVector a;
				SIMD_LOAD_A_FLOAT3_AS_0ZYX( vA, a ); 
				  
				SIMDVector result;
				SIMD_LENGTH( a, result );

				Vec3 vResult;				
				SIMD_STORE_A_FLOAT3( vResult, result ); 
				
				UTEST( vResult[0] == sqrtf(50.f) );
			}

			// --- 7. --------------------------------------------------------------------
			LOG("-- Normalize( ) \n");
			{	 		
				{
					Vec3 vA ( 3, 0, 0 );
					SIMDVector a;
					SIMD_LOAD_A_FLOAT3_AS_0ZYX( vA, a ); 
					  
					SIMDVector result;
					SIMD_NORMALIZE( a, result );

					Vec3 vResult;				
					SIMD_STORE_A_FLOAT3( vResult, result ); 
					
					UTEST( vResult == Vec3(1,0,0) );
				}
				{
					Vec3 vA ( 0, 15, 0 );
					SIMDVector a;
					SIMD_LOAD_A_FLOAT3_AS_0ZYX( vA, a ); 
					  
					SIMDVector result;
					SIMD_NORMALIZE( a, result );

					Vec3 vResult;				
					SIMD_STORE_A_FLOAT3( vResult, result ); 
					
					UTEST( vResult == Vec3(0,1,0) );
				}
				{
					Vec3 vA ( 0, 0, 87 );
					SIMDVector a;
					SIMD_LOAD_A_FLOAT3_AS_0ZYX( vA, a ); 
					  
					SIMDVector result;
					SIMD_NORMALIZE( a, result );

					Vec3 vResult;				
					SIMD_STORE_A_FLOAT3( vResult, result ); 
					
					UTEST( vResult == Vec3(0,0,1) );
				}
			}
	 		 // --- 8. --------------------------------------------------------------------
			LOG("-- LERP( ) \n");
			{	 	
				{
					Vec3 vA( 0, 1, 2 );
					Vec3 vB( 3, 4, 5 );
					float fT = 0.0f;
					SIMDVector a, b, t;
					SIMD_LOAD_A_FLOAT3( vA, a );
					SIMD_LOAD_A_FLOAT3( vB, b );
					SIMD_LOAD_A_FLOAT_AS_XXXX( fT, t );
					  
					SIMDVector result;
					SIMD_LERP( a, b, t, result );

					Vec3 vResult;				
					SIMD_STORE_A_FLOAT3( vResult, result ); 
					
					UTEST( vResult == Vec3( 0, 1, 2 ) );
				}
				{
					Vec3 vA( 0, 1, 2 );
					Vec3 vB( 3, 4, 5 );
					float fT = 1.0f;
					SIMDVector a, b, t;
					SIMD_LOAD_A_FLOAT3( vA, a );
					SIMD_LOAD_A_FLOAT3( vB, b );
					SIMD_LOAD_A_FLOAT_AS_XXXX( fT, t );
					  
					SIMDVector result;
					SIMD_LERP( a, b, t, result );

					Vec3 vResult;				
					SIMD_STORE_A_FLOAT3( vResult, result ); 
					
					UTEST( vResult == Vec3( 3, 4, 5 ) );
				}
				{
					Vec3 vA( 0, 1, 2 );
					Vec3 vB( 3, 4, 5 );
					float fT = 0.5f;
					SIMDVector a, b, t;
					SIMD_LOAD_A_FLOAT3( vA, a );
					SIMD_LOAD_A_FLOAT3( vB, b );
					SIMD_LOAD_A_FLOAT_AS_XXXX( fT, t );
					  
					SIMDVector result;
					SIMD_LERP( a, b, t, result );

					Vec3 vResult;				
					SIMD_STORE_A_FLOAT3( vResult, result ); 
					
					UTEST( vResult == Vec3( 1.5f, 2.5f, 3.5f ) );
				}
			}		
			// --- 9. --------------------------------------------------------------------
			LOG("-- Matrix - Load( ) \n");
			{	 		
				float mA[12] = {0,1,2,3,4,5,6,7,8,9,10,11};

				SIMDMatrix34 m;
				SIMD_MATRIX34_LOAD_A_12FLOAT( mA[0], mA[1], mA[2], mA[3], mA[4], mA[5], mA[6], mA[7], mA[8], mA[9], mA[10], mA[11], &m ); 
				  
				ALIGN_16( float mResult[12] );
				SIMD_STORE_A_FLOAT4( &mResult[0], m.v[0] ); 
				SIMD_STORE_A_FLOAT4( &mResult[4], m.v[1] ); 
				SIMD_STORE_A_FLOAT4( &mResult[8], m.v[2] ); 
				 
				UTEST( mResult[0] == mA[0] );					
				UTEST( mResult[1] == mA[1] );					
				UTEST( mResult[2] == mA[2] );					
				UTEST( mResult[3] == mA[3] );
				UTEST( mResult[4] == mA[4] );					
				UTEST( mResult[5] == mA[5] );					
				UTEST( mResult[6] == mA[6] );					
				UTEST( mResult[7] == mA[7] );
				UTEST( mResult[8] == mA[8] );					
				UTEST( mResult[9] == mA[9] );					
				UTEST( mResult[10] == mA[10] );					
				UTEST( mResult[11] == mA[11] );
			} 
			// --- 10. --------------------------------------------------------------------
			LOG("-- Matrix - Add( ) \n");
			{	 		
				float mA[12] = {0,1,2,3,4,5,6,7,8,9,10,11};
				float mB[12] = {0,1,2,3,4,5,6,7,8,9,10,11};

				SIMDMatrix34 m,n;
				SIMD_MATRIX34_LOAD_A_12FLOAT( mA[0], mA[1], mA[2], mA[3], mA[4], mA[5], mA[6], mA[7], mA[8], mA[9], mA[10], mA[11], &m ); 
				SIMD_MATRIX34_LOAD_A_12FLOAT( mB[0], mB[1], mB[2], mB[3], mB[4], mB[5], mB[6], mB[7], mB[8], mB[9], mB[10], mB[11], &n ); 
				
				SIMDMatrix34 r;
				SIMD_MATRIX34_ADD(&m,&n,&r);

				ALIGN_16( float mResult[12] );
				SIMD_STORE_A_FLOAT4( &mResult[0], r.v[0] ); 
				SIMD_STORE_A_FLOAT4( &mResult[4], r.v[1] ); 
				SIMD_STORE_A_FLOAT4( &mResult[8], r.v[2] ); 
				 
				UTEST( mResult[0] == mA[0] + mB[0] );
				UTEST( mResult[1] == mA[1] + mB[1] );					
				UTEST( mResult[2] == mA[2] + mB[2] );					
				UTEST( mResult[3] == mA[3] + mB[3] );
				UTEST( mResult[4] == mA[4] + mB[4] );					
				UTEST( mResult[5] == mA[5] + mB[5] );					
				UTEST( mResult[6] == mA[6] + mB[6] );					
				UTEST( mResult[7] == mA[7] + mB[7] );
				UTEST( mResult[8] == mA[8] + mB[8] );					
				UTEST( mResult[9] == mA[9] + mB[9] );					
				UTEST( mResult[10] == mA[10] + mB[10] );					
				UTEST( mResult[11] == mA[11] + mB[11] );
			} 

			// --- 11. --------------------------------------------------------------------
			LOG("-- Matrix - Mult Scalar( ) \n");
			{	 		
				float mA[12] = {0,1,2,3,4,5,6,7,8,9,10,11}; 

				SIMDMatrix34 m;
				SIMD_MATRIX34_LOAD_A_12FLOAT( mA[0], mA[1], mA[2], mA[3], mA[4], mA[5], mA[6], mA[7], mA[8], mA[9], mA[10], mA[11], &m ); 

				float fV = 2.f;
				SIMDVector v;
				SIMD_LOAD_A_FLOAT_AS_XXXX( fV, v );
				
				SIMDMatrix34 r;
				SIMD_MATRIX34_MULT_SCALAR(&m,v,&r);

				ALIGN_16( float mResult[12]);
				SIMD_STORE_A_FLOAT4( &mResult[0], r.v[0] ); 
				SIMD_STORE_A_FLOAT4( &mResult[4], r.v[1] ); 
				SIMD_STORE_A_FLOAT4( &mResult[8], r.v[2] ); 
				 
				UTEST( mResult[0] == mA[0] * fV );
				UTEST( mResult[1] == mA[1] * fV );					
				UTEST( mResult[2] == mA[2] * fV );					
				UTEST( mResult[3] == mA[3] * fV );
				UTEST( mResult[4] == mA[4] * fV );					
				UTEST( mResult[5] == mA[5] * fV );					
				UTEST( mResult[6] == mA[6] * fV );					
				UTEST( mResult[7] == mA[7] * fV );
				UTEST( mResult[8] == mA[8] * fV );					
				UTEST( mResult[9] == mA[9] * fV );					
				UTEST( mResult[10] == mA[10] * fV );					
				UTEST( mResult[11] == mA[11] * fV );
			} 

			// --- 12. --------------------------------------------------------------------
			LOG("-- Matrix - Mult Vector( ) \n");
			{	 		
				float mA[12] = {0,1,2,3,4,5,6,7,8,9,10,11}; 

				SIMDMatrix34 m;
				SIMD_MATRIX34_LOAD_A_12FLOAT( mA[0], mA[1], mA[2], mA[3], mA[4], mA[5], mA[6], mA[7], mA[8], mA[9], mA[10], mA[11], &m ); 
 				
				float vA[4] = { 0, 1, 2, 3 };
				SIMDVector v;
				SIMD_LOAD_A_FLOAT4( vA, v ); 

				SIMDVector r;
				SIMD_MATRIX34_MULT_VEC4(&m,v,r); 

				ALIGN_16( float vResult[3]);
				SIMD_STORE_A_FLOAT3( vResult, r ); 
				 
				UTEST( vResult[0] == 14 );
				UTEST( vResult[1] == 38 );
				UTEST( vResult[2] == 62 );
			} 

			return UTEST_RESULT; 
		}
	}
}
namespace Melting
{
	namespace Math
	{ 	

#ifndef __FINALROM
		void Quaternion::Print( void ) const
		{
			LOG("Real %f; Imaginary ", mfRealPart ); mvImaginaryPart.Print(); 
		}

		bool Quaternion::UnitTest( void )
		{
			LOG("::::::::: CLASS : Quaternion ::::::::::\n");

			UTEST_START;
 
	
			// --- 0. --------------------------------------------------------------------
			LOG("-- AngleAxis( Radian, Vec3 )\n");
			{	
				Quaternion q( Math::PI / 2, Vec3( 1, 0, 0 )  ); 

				UTEST( q == Quaternion ( 0.7071f, 0.7071f, 0.0f, 0.0f ) );
			}

			/*				
			// --- 1. --------------------------------------------------------------------
			LOG("-- FromMatrix( const Matrix33& )   \n");
			{	
				{
					Matrix33 mtxRotation;
					mtxRotation


					Matrix33 checkRot;

					double A    = cos( DAETools::DegreeToRadian( fAngleX ));
					double B    = sin( DAETools::DegreeToRadian( fAngleX ));
					double C    = cos( DAETools::DegreeToRadian( fAngleY ));
					double D    = sin( DAETools::DegreeToRadian( fAngleY ));
					double E    = cos( DAETools::DegreeToRadian( fAngleZ ));
					double F    = sin( DAETools::DegreeToRadian( fAngleZ ));
					double AD	=   A * D;
					double BD   =   B * D;

					checkRot.ma[0]  =   C * E;
					checkRot.ma[1]  =  -C * F;
					checkRot.ma[2]  =   D;

					checkRot.ma[3]  =  BD * E + A * F;
					checkRot.ma[4]  = -BD * F + A * E;
					checkRot.ma[5]  =  -B * C;
					checkRot.ma[6]  = -AD * E + B * F;
					checkRot.ma[7]  =  AD * F + B * E;
					checkRot.ma[8] =   A * C;
					
					UTEST( vResult == vOriginal ); 					
				} 
			}
			
			// --- 2. --------------------------------------------------------------------
			LOG("-- Add(  ) \n");
			{		 
			}
			
			*/ 
			
				
	 		return UTEST_RESULT; 
		}
#endif	 
	} // namespace Math
} // namespace Melting



namespace Melting
{	
	namespace Math
	{ 	 
#ifndef __FINALROM
 		//------------------------------------------------------------------------------ 
		// Print in the console
		//------------------------------------------------------------------------------
		void AABBox::Print( void ) const
		{
			LOG("AABBox:: [Min ="); mvMin.Print();
			LOG("; Max="); 			mvMax.Print();
			LOG("]");
		}
		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		bool AABBox::UnitTest( void )
		{
			LOG("::::::::: CLASS : AABBox ::::::::::\n");
	
			UTEST_START;
	
			// --- 0. --------------------------------------------------------------------
			LOG("-- SetupFromMinPositionAndVolume( Vec3, Vec3 )\n");
			{	
				AABBox box;
				
				box.SetupFromMinPositionAndVolume( Vec3(REAL(0)), Vec3(REAL(1)) );								
				
				UTEST( box.mvMin == Vec3(REAL(0)) && box.mvMax == Vec3(REAL(1)) );				
			}
			// --- 1. --------------------------------------------------------------------
			LOG("-- SetupFromCenterPositionAndHalfVolume( Vec3, Vec3 )\n");
			{	
				AABBox box;
				
				box.SetupFromCenterPositionAndHalfVolume( Vec3(REAL(0)), Vec3(REAL(1)) );								
				
				UTEST( box.mvMin == Vec3(REAL(-1)) && box.mvMax == Vec3(REAL(1)) );				
			}
			
			// --- 2. --------------------------------------------------------------------
			LOG("-- ComputeCenter( out Vec3 )\n");
			{	
				AABBox box(Vec3( REAL( 0 ) ), Vec3( REAL( 10 ) ));
				
				Vec3 vPoint;
				box.ComputeCenter(vPoint);
				
				UTEST( vPoint == Vec3( REAL(5) ) );				
			}
			// --- 3. --------------------------------------------------------------------
			LOG("-- IsInside( Vec3 )\n");
			{	
				AABBox box(Vec3( REAL(-1) ), Vec3( REAL(1) ));
				
				Vec3 vPoint ( REAL( 0.f ) );
				
				
				UTEST( box.IsInside( vPoint ) );				
			}
			 
			// --- 4. --------------------------------------------------------------------
			LOG("-- Union( ... )\n");
			{	
				AABBox boxA 		( Vec3( REAL(-1) ), Vec3( REAL(1) ));
				AABBox boxB 		( Vec3( REAL(0) ), Vec3( REAL(10) ));
				AABBox boxResult;
				
				AABBox::Union( boxA, boxA, boxResult );
				UTEST( boxResult == AABBox( Vec3( REAL(-1) ), Vec3( REAL(1))));
				
				AABBox::Union( boxA, boxB, boxResult );
				UTEST( boxResult == AABBox( Vec3( REAL(-1) ), Vec3( REAL(10))));
				
				AABBox::Union( boxA, Vec3( REAL(2)), boxResult );
				UTEST( boxResult == AABBox( Vec3( REAL(-1) ), Vec3( REAL(2))));
			} 			
			
			
			
	 		return UTEST_RESULT; 
		}
#endif
 
	}
}

namespace Melting
{
	namespace Math
	{ 	 
	
#ifndef __FINALROM
		 			
		// Print infos to the console
		void BoundingSphere::Print( void ) const
		{
			LOG("Sphere [ Center = ");		mvCenter.Print();			
#if 0
			LOG("; Radius = ");				mfxRadius.Print();
#endif
			LOG("]"); 
			
//				LOG("Sphere::Radius² = \n");mfxSquaredRadius.Print();
		}   

		bool BoundingSphere::UnitTest( void )
		{
			LOG("::::::::: CLASS : BoundingSphere ::::::::::\n");
	
			UTEST_START;
	
			// --- 0. Operator == (required for most tests) --------------------------------------------------------------------
			LOG("-- Operator == ...\n");
			{	
				UTEST( BoundingSphere ( Vec3( REAL(0)), REAL(2)) == BoundingSphere ( Vec3( REAL(0)), REAL(2)));				
				UTEST( BoundingSphere ( Vec3( REAL(1)), REAL(2)) == BoundingSphere ( Vec3( REAL(1)), REAL(2)));				
				UTEST( BoundingSphere ( Vec3( REAL(0), REAL(1) , REAL(2)), REAL(2)) == BoundingSphere ( Vec3( REAL(0), REAL(1) , REAL(2)), REAL(2)));				
			}
	
#if 0		
			// --- 1. SetupFromAABBox --------------------------------------------------------------------
			LOG("-- SetupFromAABBox...\n");
			{ 
				AABBox box( Vec3( REAL( 0.f ) ), Vec3( REAL( 10 ) ));
								
				BoundingSphere bSphere;
				
				Vec3 vCenter ( REAL(5) );
				
				bSphere.SetupFromAABBoxInside( box );
				UTEST( BoundingSphere ( vCenter, REAL(75 /* = 5² par dimension */).SquaredRoot() ) == bSphere );	
				
				bSphere.SetupFromAABBoxOutside( box );
				UTEST( BoundingSphere ( Vec3( REAL(5)), REAL(5) ) == bSphere );
			}
			
#endif
				
			// --- 2. IsInside --------------------------------------------------------------------
			LOG("-- IsInside...\n");
			{ 
				BoundingSphere 	boundingSphereA ( Vec3( REAL(-1), REAL(0), REAL(0)), REAL(3));
			
				
				UTEST( boundingSphereA.IsInside		( Vec3 ( REAL(0), REAL(0), REAL(0) ))) ;
				UTEST( boundingSphereA.IsInside		( Vec3( REAL(0), REAL(1), REAL(0) ) )) ;
				UTEST( ! boundingSphereA.IsInside	( Vec3( REAL(10), REAL(0), REAL(0) ) )) ;
				UTEST( ! boundingSphereA.IsInside	( Vec3( REAL(0), REAL(), REAL(5) ) )) ;
			}
			
			
			// --- 3. Union --------------------------------------------------------------------
			LOG("-- Union ...\n");
			{
				BoundingSphere boundingSphereA ( Vec3( REAL(-1), REAL(0), REAL(0)), REAL(1));
				BoundingSphere boundingSphereB ( Vec3( REAL(1), REAL(0), REAL(0)), REAL(1)); 
				BoundingSphere boundingSphereUnionResult;	
				BoundingSphere::Union( boundingSphereA , boundingSphereB, boundingSphereUnionResult );
				
				UTEST( boundingSphereUnionResult == BoundingSphere ( Vec3( REAL(0)), REAL(2)));				
			}
			
			{
				BoundingSphere boundingSphereA ( Vec3( REAL(0), REAL(-1), REAL(0)), REAL(1));
				BoundingSphere boundingSphereB ( Vec3( REAL(0), REAL(1), REAL(0)), REAL(1)); 
				BoundingSphere boundingSphereUnionResult;	
				BoundingSphere::Union( boundingSphereA , boundingSphereB, boundingSphereUnionResult );
				
				UTEST( boundingSphereUnionResult == BoundingSphere ( Vec3( REAL(0)), REAL(2)));				
			}
			{
				BoundingSphere boundingSphereA ( Vec3( REAL(0), REAL(0), REAL(-1)), REAL(1));
				BoundingSphere boundingSphereB ( Vec3( REAL(0), REAL(0), REAL(1)), REAL(1)); 
				BoundingSphere boundingSphereUnionResult;	
				BoundingSphere::Union( boundingSphereA , boundingSphereB, boundingSphereUnionResult );
				
				UTEST( boundingSphereUnionResult == BoundingSphere ( Vec3( REAL(0)), REAL(2)));				
			}
			
			{
				BoundingSphere boundingSphereA ( Vec3( REAL(-0.1f), REAL(0), REAL(0)), REAL(0.1f));
				BoundingSphere boundingSphereB ( Vec3( REAL(0.1f), REAL(0), REAL(0)), REAL(0.1f)); 
				BoundingSphere boundingSphereUnionResult;	
				BoundingSphere::Union( boundingSphereA , boundingSphereB, boundingSphereUnionResult );
				
				UTEST( boundingSphereUnionResult == BoundingSphere ( Vec3( REAL(0)), REAL(0.2f)));				
			}
	
			{
				BoundingSphere boundingSphereA ( Vec3( REAL(0), REAL(-0.1f), REAL(0)), REAL(0.1f));
				BoundingSphere boundingSphereB ( Vec3( REAL(0), REAL(0.1f), REAL(0)), REAL(0.1f)); 
				BoundingSphere boundingSphereUnionResult;	
				BoundingSphere::Union( boundingSphereA , boundingSphereB, boundingSphereUnionResult );
				
				UTEST( boundingSphereUnionResult == BoundingSphere ( Vec3( REAL(0)), REAL(0.2f)));				
			}
	
	
			{
				BoundingSphere boundingSphereA ( Vec3( REAL(0), REAL(0), REAL(-0.1f) ), REAL(0.1f));
				BoundingSphere boundingSphereB ( Vec3( REAL(0), REAL(0), REAL(0.1f) ), REAL(0.1f)); 
				BoundingSphere boundingSphereUnionResult;	
				BoundingSphere::Union( boundingSphereA , boundingSphereB, boundingSphereUnionResult );
				
				UTEST( boundingSphereUnionResult == BoundingSphere ( Vec3( REAL(0)), REAL(0.2f)));				
			}
			
			{
				BoundingSphere boundingSphereA ( Vec3( REAL(1), REAL(0), REAL(0)), REAL(1));
				BoundingSphere boundingSphereB ( Vec3( REAL(9), REAL(0), REAL(0)), REAL(3)); 
				BoundingSphere boundingSphereUnionResult;	
				BoundingSphere::Union( boundingSphereA , boundingSphereB, boundingSphereUnionResult );
				
				UTEST( boundingSphereUnionResult == BoundingSphere ( Vec3( REAL(6), REAL(0), REAL(0)), REAL(6)));				
			}
			
			{
				BoundingSphere boundingSphereA ( Vec3( REAL(0), REAL(1), REAL(0)), REAL(1));
				BoundingSphere boundingSphereB ( Vec3( REAL(0), REAL(9), REAL(0)), REAL(3)); 
				BoundingSphere boundingSphereUnionResult;	
				BoundingSphere::Union( boundingSphereA , boundingSphereB, boundingSphereUnionResult );
				
				UTEST( boundingSphereUnionResult == BoundingSphere ( Vec3( REAL(0), REAL(6), REAL(0)), REAL(6)));				
			}
			
			{
				BoundingSphere boundingSphereA ( Vec3( REAL(0), REAL(0), REAL(1) ), REAL(1));
				BoundingSphere boundingSphereB ( Vec3( REAL(0), REAL(0), REAL(9) ), REAL(3)); 
				BoundingSphere boundingSphereUnionResult;	
				BoundingSphere::Union( boundingSphereA , boundingSphereB, boundingSphereUnionResult );
				
				UTEST( boundingSphereUnionResult == BoundingSphere ( Vec3( REAL(0), REAL(0), REAL(6)), REAL(6)));				
			}
			
			{ // different case : B is included in A !
			
				BoundingSphere boundingSphereA ( Vec3( REAL(0) ), REAL(2));
				BoundingSphere boundingSphereB ( Vec3( REAL(0), REAL(0), REAL(1) ), REAL(0.5f)); 
				BoundingSphere boundingSphereUnionResult;	
				BoundingSphere::Union( boundingSphereA , boundingSphereB, boundingSphereUnionResult );
				
				UTEST( boundingSphereUnionResult == BoundingSphere ( Vec3( REAL(0) ), REAL(2) ));				
			}
			
	 		return UTEST_RESULT;
		} 
#endif
	}
}
namespace Melting
{
	namespace Math
	{ 	 
 		
#ifndef __FINALROM		
			//------------------------------------------------------------------------------ 
			// Print in the console
			//------------------------------------------------------------------------------
			void SRT::Print( void ) const
			{
				LOG("SRT:: [Scale="); 	mvScale.Print();
				LOG("; Rotate="); 		mquatRotation.Print();
				LOG("; Translate=");	mvTranslate.Print(); 
				LOG("]");
			}
			
			//------------------------------------------------------------------------------
			//------------------------------------------------------------------------------
			bool SRT::UnitTest( void )
			{ 
				LOG("::::::::: CLASS : SRT ::::::::::\n");
		
				UTEST_START; 

				// --- 0. --------------------------------------------------------------------
				LOG("void SRT::BuildIdentity( SRT& _rOutIdentity)\n");
				{
					SRT		ID( Vec3::ONE, Quaternion(1,0,0,0), Vec3::ZERO );

					SRT		result;
					SRT::BuildIdentity( result );

					UTEST( result == ID );			
				}

				// --- 1. --------------------------------------------------------------------
				LOG("void SRT::Transform( const SRT& _rLocal, const SRT& _rParent, SRT& _rOutGlobal )\n");
				{
					{
						LOG(". Parent is Identity\n");

						SRT		parent;
						SRT::BuildIdentity( parent );

						SRT local(Vec3(5,1,1), Quaternion(2,Vec3(1,0,0)),Vec3(3,2,1));

						SRT result;

						SRT::Transform( local, parent, result );

						UTEST( result == local );
					}
					{
						LOG(". Parent rotation is 90° around UNIT_Y, Local is vector UNIT_X, so we should have Vec3(0,0,-1)\n");

						SRT	parent( Vec3::ONE, Quaternion( Math::DegreeToRadian( 90 ), Vec3(0,1,0)), Vec3::ZERO);

						SRT local( Vec3::ONE, Quaternion(0,Vec3(1,0,0)),Vec3(1,0,0) );

						SRT result;

						SRT::Transform( local, parent, result );

						UTEST( result.mvTranslate == Vec3(0,0,-1) );
					}
				}
 
				
		 		return UTEST_RESULT; 
			}
#endif	
	}
}
namespace Melting
{
	namespace Math
	{ 	 
#ifndef __FINALROM
 		//------------------------------------------------------------------------------ 
		// Print in the console
		//------------------------------------------------------------------------------
		void Cone::Print( void ) const
		{
#if 0
			LOG("Cone::Vertex ="); 			mvVertex.Print();			LOG("\n");
			LOG("Cone::Direction ="); 		mvDirection.Print();		LOG("\n");
			LOG("Cone::Length ="); 			mfxLength.Print();			LOG("\n");
			LOG("Cone::HalfAngleCos ="); 	mfxHalfAngleCos.Print();	LOG("\n");
			LOG("Cone::HalfAngleSin ="); 	mfxHalfAngleSin.Print(); 	LOG("\n");
#endif
		}
		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		bool Cone::UnitTest( void )
		{
			LOG("::::::::: CLASS : Cone ::::::::::\n");
	
			UTEST_START;
	
			// --- 0. --------------------------------------------------------------------
			LOG("-- SetupFromTwoPointsAndAngle( ... )\n");
			{	
				Cone cone;
				cone.SetupFromTwoPointsAndAngle( Vec3(REAL( 0.f )),  Vec3(REAL( 0.f ), REAL( 0.f ), REAL( -100 ) ), REAL( 45 ) ); 
			
				UTEST( cone == Cone( 	Vec3( REAL( 0.f ) ) , Vec3(REAL( 0.f ), REAL( 0.f ), REAL( -1 )) , REAL(100), REAL(45) ) );
			}
			
			// --- 0. --------------------------------------------------------------------
			LOG("-- IsInside( Vec3 )\n");
			{	
				Cone cone;
				cone.SetupFromTwoPointsAndAngle( Vec3(REAL( 0.f )),  Vec3(REAL( 0.f ), REAL( 0.f ), REAL( -100 ) ), REAL( 45 ) );
				
				
				UTEST( Cone::IsInside( Vec3 ( REAL( 0.f ), REAL( 0.f ), REAL( -40 ) ) , cone ) );				
				UTEST( ! Cone::IsInside( Vec3 ( REAL( 0.f ), REAL( 0.f ), REAL( 40 ) ) , cone ) );				
			}
			
			
	 		return UTEST_RESULT; 
		}
#endif
	}
}

namespace Melting
{
	namespace Math
	{	
	
		
#ifndef __FINALROM
		//------------------------------------------------------------------------------ 
		// Print in the console
		//------------------------------------------------------------------------------
		void Plane::Print( void ) const
		{
			LOG("Plane:: [Normal ="); 	mvNormal.Print();
#if 0
			LOG(";Dist ="); 	mfxDist.Print();
#endif
			LOG("]");
		}
		
		//------------------------------------------------------------------------------
		//------------------------------------------------------------------------------
		bool Plane::UnitTest( void )
		{
			LOG("::::::::: CLASS : Plane ::::::::::\n");
			
			UTEST_START; 

			// --- 0. --------------------------------------------------------------------
			LOG("-- SetupFromTriangle(...) && SetupFromNormalAndPoint(...) \n");
			{
				Plane plane;
				plane.SetupFromTriangle( Vec3(REAL(1), REAL(0), REAL(0)), Vec3(REAL(0), REAL(0), REAL(0)), Vec3(REAL(0), REAL(0), REAL(1)) );				
				UTEST( Plane( VectorUp(), REAL( 0.f ) ) == plane );
				
				plane.SetupFromNormalAndPoint( Vec3(REAL(1), REAL(0), REAL(0)), Vec3( REAL( 0.f ) ) );
				UTEST( Plane( VectorRight(), REAL( 0.f ) ) == plane );
			}
			
			
			// --- 1. --------------------------------------------------------------------
			LOG("-- DistanceFrom \n");
			{
				REAL fxDist (10);
				
				Plane plane ( Vec3( REAL( 0.f ), REAL( 1.f ), REAL( 0.f ) ), REAL( 0.f ) );
				
				Vec3 vPoint ( REAL( 0.f ), fxDist, REAL( 0.f ) );
				
				
				UTEST( fxDist == plane.DistanceFrom( vPoint ) );
			}
			
			// --- 2. --------------------------------------------------------------------
			LOG("-- Project \n");
			{
				Plane plane ( VectorUp(), REAL(1) );
				
				Vec3 vPoint ( REAL( 0.f ), REAL(10), REAL( 0.f ) );
				Vec3 vProjected;
				
				plane.Project( vPoint, vProjected );
				
				UTEST( Vec3(REAL(0), REAL(1), REAL(0)) == vProjected );
			}
			
			// --- 3. --------------------------------------------------------------------
			LOG("-- IsAbove \n");
			{
				Plane plane( VectorUp(), REAL(1) );
				
				UTEST( plane.IsAbove( Vec3( REAL(2) ) ) );
				
				UTEST( ! plane.IsAbove( Vec3( REAL(0) ) ) );
			}
			
			
	 		return UTEST_RESULT; 
		}
		
#endif		 
	} // namespace Math
		
} // namespace Melting
 
namespace Melting
{
	namespace Math
	{ 	    
			
#ifndef __FINALROM		
			//------------------------------------------------------------------------------ 
			// Print in the console
			//------------------------------------------------------------------------------
			void Ray::Print( void ) const
			{
				LOG("Ray:: [Origin="); 			mvOrigin.Print();
				LOG("; Direction="); 			mvDirection.Print();
				LOG("; InverseDirection =");	mvInverseDirection.Print(); 
				LOG("]");
			}
			
			//------------------------------------------------------------------------------
			//------------------------------------------------------------------------------
			bool Ray::UnitTest( void )
			{ 
				LOG("::::::::: CLASS : Ray ::::::::::\n");
		
				UTEST_START;
		
				/*
				// --- 0. --------------------------------------------------------------------
				LOG("-- SetupFromSegment \n");
				{	
					Ray 	ray;
					ray.SetupFromSegment( Vec3::ZERO, Vec3( REAL(0), REAL(10), REAL(0) ) )
					
					UTEST( ray = Ray(Vec3::ZERO, Vec3(REAL(0),Vec3(REAL(1),Vec3(REAL(0)));
				}
				*/
				// --- 1. --------------------------------------------------------------------
				LOG("-- ComputePoint(  ) \n");
				{	 
					Ray 	ray( Vec3( REAL(0), REAL(0), REAL(10) ), VectorBack() );
					
					Vec3 vResult;
					
					ray.ComputePoint( REAL(10) , vResult);
					UTEST( vResult == Vec3::ZERO );					
				}				
				
				// --- 2. --------------------------------------------------------------------
				LOG("-- ComputeDistance(  ) \n");
				{	 
					Ray 	ray( Vec3( REAL(0), REAL(0), REAL(10) ), VectorBack() );
					Plane 	plane( VectorFront(), REAL(1) ) ;
					
					UTEST( ray.ComputeDistance( plane ) == REAL(9) );					
				}
				
				// --- 3. --------------------------------------------------------------------
				LOG("-- ComputeIntersection(  ) \n");
				{	 
					Ray 	ray( Vec3( REAL(0), REAL(0), REAL(10) ), VectorBack() );
					Plane 	plane( VectorFront(), REAL(5) ) ;
					
					Vec3 vResult;
					ray.ComputeIntersection( plane , vResult );
					
					UTEST( vResult == Vec3( REAL(0), REAL(0), REAL(5) ) );					
				}
				// --- 4. --------------------------------------------------------------------
				LOG("-- DoesIntersect( Plane ) \n");
				{	 				
					Ray 	ray( Vec3( REAL(0), REAL(0), REAL(10) ), VectorBack() );	
					Plane 	plane( Vec3( REAL(1) ).Normalize(), REAL( 0.f ) );					
					
					UTEST( ray.DoesIntersect( plane ) );
					
					plane = Plane (VectorUp(), REAL( 0.f ) );
					UTEST( ! ray.DoesIntersect( plane ) );
				}
				
#if 0
				// --- 5. --------------------------------------------------------------------
				LOG("-- DoesIntersect( AABBox ) \n");
				{
					
#define 	TestIntersectRayAABBox(boxCenter,rayOrigin)										\
					box = AABBox ( boxCenter - vBoxHalfSize, boxCenter + vBoxHalfSize );	\
					UTEST( Ray ( rayOrigin, (boxCenter-rayOrigin).Normalize()).DoesIntersect( box ) );
					
					Ray 	ray;
					Vec3 	vBoxCenter	;
					Vec3	vBoxHalfSize	( REAL(3) );					
					
					
					AABBox 	box( vBoxCenter - vBoxHalfSize, vBoxCenter + vBoxHalfSize );					
					
					LOG("RAY FROM OUTSIDE THE BOX!!!\n");
					LOG("BoxCenter = (0,0,0)\n");
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(0), REAL(10) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(0), REAL(-10) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(-10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(0)), 	Vec3( REAL(10), REAL(0), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(0)), 	Vec3( REAL(-10), REAL(0), REAL(0) )); 	
					 
					
					LOG("BoxCenter = (5,0,0)\n");
					TestIntersectRayAABBox(Vec3( REAL(5), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(0), REAL(10) ));
					TestIntersectRayAABBox(Vec3( REAL(5), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(0), REAL(-10) ));
					TestIntersectRayAABBox(Vec3( REAL(5), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(5), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(-10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(5), REAL(0), REAL(0)), 	Vec3( REAL(10), REAL(0), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(5), REAL(0), REAL(0)), 	Vec3( REAL(-10), REAL(0), REAL(0) ));  
					
					LOG("BoxCenter = (-5,0,0)\n");
					TestIntersectRayAABBox(Vec3( REAL(-5), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(0), REAL(10) ));
					TestIntersectRayAABBox(Vec3( REAL(-5), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(0), REAL(-10) ));
					TestIntersectRayAABBox(Vec3( REAL(-5), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(-5), REAL(0), REAL(0)), 	Vec3( REAL(0), REAL(-10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(-5), REAL(0), REAL(0)), 	Vec3( REAL(10), REAL(0), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(-5), REAL(0), REAL(0)), 	Vec3( REAL(-10), REAL(0), REAL(0) )); 						
				  
															
					
					LOG("BoxCenter = (0,5,0)\n");
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(5), REAL(0)), 	Vec3( REAL(0), REAL(0), REAL(10) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(5), REAL(0)), 	Vec3( REAL(0), REAL(0), REAL(-10) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(5), REAL(0)), 	Vec3( REAL(0), REAL(10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(5), REAL(0)), 	Vec3( REAL(0), REAL(-10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(5), REAL(0)), 	Vec3( REAL(10), REAL(0), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(5), REAL(0)), 	Vec3( REAL(-10), REAL(0), REAL(0) )); 
					
					LOG("BoxCenter = (0,-5,0)\n");
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(-5), REAL(0)), 	Vec3( REAL(0), REAL(0), REAL(10) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(-5), REAL(0)), 	Vec3( REAL(0), REAL(0), REAL(-10) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(-5), REAL(0)), 	Vec3( REAL(0), REAL(10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(-5), REAL(0)), 	Vec3( REAL(0), REAL(-10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(-5), REAL(0)), 	Vec3( REAL(10), REAL(0), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(-5), REAL(0)), 	Vec3( REAL(-10), REAL(0), REAL(0) )); 
					  
											
				 
					LOG("BoxCenter = (0,0,5)\n");
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(5)), 	Vec3( REAL(0), REAL(0), REAL(10) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(5)), 	Vec3( REAL(0), REAL(0), REAL(-10) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(5)), 	Vec3( REAL(0), REAL(10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(5)), 	Vec3( REAL(0), REAL(-10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(5)), 	Vec3( REAL(10), REAL(0), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(5)), 	Vec3( REAL(-10), REAL(0), REAL(0) ));
				 
					LOG("BoxCenter = (0,0,-5)\n");
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(-5)), 	Vec3( REAL(0), REAL(0), REAL(10) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(-5)), 	Vec3( REAL(0), REAL(0), REAL(-10) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(-5)), 	Vec3( REAL(0), REAL(10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(-5)), 	Vec3( REAL(0), REAL(-10), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(-5)), 	Vec3( REAL(10), REAL(0), REAL(0) ));
					TestIntersectRayAABBox(Vec3( REAL(0), REAL(0), REAL(-5)), 	Vec3( REAL(-10), REAL(0), REAL(0) )); 	 					
					// RAY FROM INSIDE THE BOX!!!
					
#undef 	TestIntersectRayAABBox	
					
				}
#endif
				
				// --- 6. --------------------------------------------------------------------
				LOG("-- DoesIntersect( BoundingSphere ) \n");
				{	 
					BoundingSphere 	sphere( Vec3::ZERO, REAL( 1 ) );
					Ray 			ray( Vec3( REAL(0), REAL(0), REAL(10) ), VectorBack() );					
					
					UTEST( ray.DoesIntersect( sphere ) );
				}		
				
		 		return UTEST_RESULT; 
			}
#endif	
	}
}
#endif // __MELTINGMATH_H__